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ABSTRACT 

We present the new code ALCAR developed to model multidimensional, multi energy- 
group neutrino transport in the context of supernovae and neutron-star mergers. The 
algorithm solves the evolution equations of the Oth- and lst-order angular moments 
of the specific intensity, supplemented by an algebraic relation for the 2nd-moment 
tensor to close the system. The scheme takes into account frame-dependent effects 
of order 0(y/c ) as well as the most important types of neutrino interactions. The 
transport scheme is significantly more efficient than a multidimensional solver of the 
Boltzmann equation, while it is more accurate and consistent than the flux-limited 
diffusion method. The finite-volume discretization of the essentially hyperbolic sys¬ 
tem of moment equations employs methods well-known from hydrodynamics. For the 
time integration of the potentially stiff moment equations we employ a scheme in 
which only the local source terms are treated implicitly, while the advection terms 
are kept explicit, thereby allowing for an efficient computational parallelization of the 
algorithm. We investigate various problem setups in one and two dimensions to verify 
the implementation and to test the quality of the algebraic closure scheme. In our 
most detailed test, we compare a fully dynamic, one-dimensional core-collapse sim¬ 
ulation with two published calculations performed with well-known Boltzmann-type 
neutrino-hydrodynamics codes and we find very satisfactory agreement. 

Key words: radiative transfer - neutrinos - hydrodynamics - supernovae - neutron 
stars 


1 INTRODUCTION 


In various astrophysical scenarios involving matter in a hot 
and dense phase, neutrino interactions take place in a way 
that the full transport problem - which consistently fol¬ 
lows the emission, propagation and absorption of neutrinos 
- needs to be taken into account to correctly describe these 
systems. A prominent example is a core-collapse supernova 
(CCSN), in which according to the present standard model 
the explosion is essentially only made possible by the en¬ 
ergy deposition due to the re-absorption of neutrinos just 


produced in the proto-neutron star (proto-NS; see Janka 
2012 Burrows]|2013 for recent reviews). Genuine neutrino- 
transport effects can also be crucial for determining the 
properties of potentially nucleosynthesis-relevant outflows 
and may even give rise to these outflows to begin with. Such 
outflows are believed to occur during a CCSN in the form 
of a neutrino-driven wind expelled from the surface of the 
proto-NS (e.g. Qian & WoosleyTl99 6). Another example is a 
massive NS formed during the merger of two NSs (Dessart 


|et al.| |2009), or similarly a black-hole (BH) torus configura¬ 
tion also produced by such a merger or by the merger of a 
NS and a BH (e.g Wanajo fe Janka|2012 Fernandez & Met- 


|zger|2013| ) . Another astrophysical scenario in which neutrino 
transport may be crucial is the launching of a gamma-ray 
burst jet, which could be powered to some degree by neu¬ 
trino pairs annihilating in the polar regions of a BH-torus 
system (e.g. Popham et al. 1999). 


Unfortunately, most multidimensional results for the 
aforementioned scenarios stem from more or less idealized 
investigations, mainly owing to the enormous computational 
requirements of a time-dependent, multidimensional treat¬ 
ment of neutrino transport. The level of simplification is typ¬ 
ically chosen to provide the optimal balance between accu¬ 
racy and computational expense, given the constraints of the 
available computational resources and the considered phys¬ 
ical effects to be captured to a sufficient degree. The most 
accurate neutrino-transport schemes follow the full spatial, 
energetic and directional dependence of the neutrino dis- 
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tribution, described by the Boltzmann equation. However, 
multidimensional applications of these “Boltzmann solvers” 
are severely constrained by their complexity and computa¬ 
tional expense and force modern simulations to employ mis¬ 
cellaneous restrictions, such as the “ray-by-ray” approach 


(e.g. Buras et al. 2006 Hanke et al. 20121, the omission 
of (a subset of) velocity-dependent terms and the decou¬ 


pling of neutrinos with different energies (e.g. Livne et al. 
2004 Qtt et al.|2008 I, or the investigation of only individual 


static background configurations of matter (Sumiyoshi et al. 


20151. A computationally cheaper alternative to a Boltz¬ 


mann solver is the “flux-limited diffusion” (FLD) method 
(see, e.g., Levermore fe Pomraning||1981 Mihalas & Miha- 


la|[l984|, in which the radiation energy density (as repre¬ 


sented by the 0th angular moment of the specific intensity, 
see Sec. [2] for explicit definitions) is evolved assuming that 
the radiation flux density (as represented by the 1st angu¬ 
lar moment) can be written as a function of the radiation 
energy density. This method is the simplest realization of a 
“truncated moment scheme”. However, FLD implies several 
physical as well as technical drawbacks. 

In this paper we present a genuinely multidimensional, 
fully energy-dependent radiation-hydrodynamics (RHD) 
scheme where both the energy density and the flux density 
of neutrinos are evolved quantities. The system of equations 
is closed by an algebraic expression for the Eddington ten¬ 
sor (defined as the normalized 2nd angular moment of the 
specific intensity) as function of the lower-order moments. 
This makes the scheme a specific realization of the class 
of “variable Eddington factor” methods and it will there¬ 
fore be denoted here as “algebraic Eddington factor” (AEF) 
method. 

In both methods, AEF and FLD, the closure (i.e. the 
corresponding angular moment which closes the set of mo¬ 
ment equations) is assumed to be given solely as a function 
of the local values of the evolved radiation moments (or ad¬ 
ditionally of the local gas properties in the case of FLD). 
Since in general the closure is a non-local function of the 
surrounding radiation sources, it is clear that a universal 
closure relation between the local angular moments, corre¬ 
sponding to an equation of state as in hydrodynamics, does 
not exist for arbitrarily shaped radiation fields. It is nev¬ 
ertheless conceivable that for not too complex geometries 
(e.g. single-source configurations like a PNS or a NS-merger 
remnant) the radiation field tends to be arranged such that 
typical relations between the lowest angular moments are 
fulfilled up to a sufficient degree. Keeping also in mind that 
not the full angular information but rather only the angular 
moments of the radiation field enter the hydrodynamics part 
of the full RHD system, the fluid may only be secondarily 
affected by the error introduced by the approximate closure. 
Certainly, the actual performance and applicability of either 
method, FLD or AEF, is problem dependent and has to be 
examined individually for a given problem setup. However, 
the main advantage of an AEF scheme (and also of an FLD 
scheme) compared to a more accurate Boltzmann solver is 
its computational efficiency, which allows to perform gen¬ 
uinely multidimensional RHD simulations and larger sets of 
model calculations with reasonable computational effort in 
the first place. 

Compared to a standard FLD solver, the AEF method 
as presented in this paper features the following advantages: 


(1) It is potentially more accurate, simply on account of the 
fact that an FLD scheme is essentially an AEF scheme in 
which a certain collection of terms is dropped (e.g. Lev- 


ermore fc Pomraning||1981| |Cernohorsky fc van den 


lorn 


1990 


Dgani & Cernohorsky 1991). (2) It is more consis¬ 


tent: A particular consequence of retaining the evolution 
equations for the 1st moments in AEF is that the conserva¬ 
tion of the total (radiation plus fluid) momentum and there¬ 
fore also of the total energy can be ensured in the case of 
momentum-exchanging interactions between neutrinos and 


the gas (Baron et al. 1989 Cernohorsky & van den Horn 


19901. Another, related advantage is that the lst-moment 


vectors can in principle point into arbitrary directions in an 
AEF scheme and are not forced to be parallel to the gradient 
of the 0th moments as in FLD. This allows, for example, to 
describe shadows behind illuminated opaque structures (see 


the test problem in Sec. 4.2.11. (3) The different mathe¬ 


matical nature of the evolution equations in AEF (hyper¬ 
bolic as opposed to parabolic in FLD) enables the applica¬ 
tion of time-explicit methods based on Riemann-solvers for 
the advection part of the system. Such methods are well- 
known from hydrodynamics. In contrast, in a time-explicit 
treatment of FLD the time step would in principle be un¬ 
bound from below, which in practice forces one to employ a 
fully time-implicit scheme. An explicit scheme compared to 
an implicit one is particularly advantageous in multidimen¬ 
sional simulations, however, because first, its computational 
requirements only scale linearly with the size of the grid, and 
second, the computational parallelization is generally more 
efficient and straightforward. 

The strategy of using an AEF scheme for radiation 
transport is not nev0 Multidimensional applications regard¬ 
ing photon transport exist in a number of realizations (see, 


e-g., 

Audit et al. 2002; 

Hayes & Norman 2003| Gonzalez 

et al. 

p007. Aubert b r 

'eyssier 2008 Vaytet et al. 2011 

Sadowski et al. 

2013; Skinner & Ostriker 2013[ |McKin- 

ney et al. 

2014 

. Until recently, however, only a few in- 


vestigations considered the AEF scheme in the context of 


neutrino transport (Schinder & Bludman 


1989 


Dgani & 


Cernohorsky| 

1991J Koerner Sz Janka 1992; |Cernohorsky 

& van Weert 

1992) , although several studies concerning 


aspects of the closure prescription (Cernohorsky & Blud 
|man|1994j |Bludman fe Cernohorsky|1995 |Smit et al.|2000| 


and the solution strategies (Smit et al. 1997 Pons et al. 


20001 elucidated its capabilities. In recent works by Shi- 
bata et al.| ( |2011| ) and |Cardall et al.| ( |2013| ) the truncated 
two-moment scheme was extended to a general relativis¬ 


tic framework. Shibata & Sekiguchi (2012) made use of 


this framework in the “grey” approximation (i.e. averag¬ 
ing over the neutrino energies) and presented axisymmet- 
ric neutrino-magnetohydrodynamics simulations of BH ac¬ 
cretion tori as models for central engines of gamma-ray 


bursts. Kuroda et al. (2012 20141 combined the relativis¬ 


tic AEF scheme with an isotropic diffusion source approx¬ 
imation (IDSA, as developed in |Liebendorfer et al.||2009| ) 
to perform three-dimensional CCSN simulations. A gen- 


1 In the literature, these schemes are sometimes simply denoted 
as “Mi schemes”, which actually only refers to the specific M\ 
closure being used to express the Eddington tensor in a truncated 


two-moment scheme (see Sec. 2.4.2 for details). 
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eral relativistic AEF scheme was implemented in the GRID 
code by O’Connor & Ott (20131, who conducted one¬ 
dimensional, energy-dependent CCSN simulations neglect¬ 
ing all velocity-dependent terms and energy-coupling inter¬ 
actions. Recently, O’Connor (2014) released an improved 
version of this scheme, including velocity-dependent terms 
and energy-coupling interactions, but still constrained to 
spherical symmetry. Another implementation of the AEF 
scheme, in special relativity and using a grey approxima¬ 
tion, was presented in |Takaha shi et al. (20131. 

In contrast to the schemes used in the aforementioned 
studies, the neutrino-transport code presented here, which 
will be designated the name “ALCAR” (Algebraic Local 
Closure Approach for Radiation), combines all of the follow¬ 
ing features: It is genuinely multidimensional, fully energy 
dependent and it takes into account all velocity-dependent 
terms up to order 0{v/c) as well as energy-coupling due 
to Doppler shift and inelastic neutrino-matter interactions. 
To validate the numerical implementation of the algorithm 
and each of its aforementioned features we examine several 
test problems. However, the test problems are not only con¬ 
ducted to check the correct numerical implementation of 
the algorithm but also serve to examine the central approx¬ 
imation of the scheme embodied in the algebraic closure. 
In particular, two tests focus on the question how the AEF 
method performs in a CCSN setup and how sensitive the 
results are to specific choices of the closure. In one test 
problem, only the neutrino field is evolved in a fixed hydro- 
dynamic background, making it possible to test the AEF 
scheme against the FLD and a Boltzmann scheme indepen¬ 
dently of the hydrodynamic treatment. In the other test 
problem, the results of a fully dynamic CCSN calculation 
are compared with results published in |Lieb endorfer et al.| 
(20051, which were obtained with the two well-established 


Boltzmann-RHD codes VERTEX-PROMETHEUS iRampp 


& Janka 20021 and AGILE-BOLTZTRAN (Liebendorfer 


et al. 2004). Our main findings from these tests are that 


the AEF scheme is sufficiently accurate to represent a com¬ 
petitive, though computationally much simpler, alternative 
to Boltzmann solvers in simulations of CCSNe. Note that 
this conclusion is in agreement with O’Connor (20141, who 
conducted a similar comparison between ID CCSN simula¬ 
tions performed with his general relativistic AEF code and 
with the two aforementioned Boltzmann codes. 

This article is structured as follows: In Sec. [2] we outline 
the steps that lead from the equation of radiative transfer to 
the moment equations that are solved in our code, and we 
motivate and present the approximations and assumptions 
contained in the AEF method. Subsequently, in Sec. [3] we 
describe the methods used to discretize the moment equa¬ 
tions in space, energy and time. In Sec. [4] we present a suite 
of test problems, in which we investigate the quality of the 
AEF approximation and study features such as the velocity 
dependence, the correct behavior in the static and dynamic 
diffusion limits, and the ability of the scheme to describe 
radiation shadowing. Finally, in Sec. [5] we summarize our 
presentation. 

The conventions regarding our notation are as follows: 
We use lower-case, italic letters to denote spatial 

tensor components, lower-case roman letters i, j, k .. . for grid 
indices and n for the time index. Moreover, we make use of 
the Einstein notation to write sums of products of tensor 


components. Symbols with hats, as for instance X, refer to 
discretized quantities, while symbols with bars, as for in¬ 
stance X, denote versions of the corresponding quantities 
integrated over the whole energy spectrum. Vectors in spa¬ 
tial and momentum space are denoted as x and p, and t 
refers to the time coordinate. The symbols c, h and fee mean 
the speed of light, the Planck constant and the Boltzmann 
constant, respectively. 


2 THE 0{V/C) EQUATIONS OF RADIATION 
HYDRODYNAMICS 

In this section we briefly define the basic quantities and 
present the RHD equations used in our code. 


2.1 The equation of radiative transfer in the 
comoving frame 

Both the equations of hydrodynamics and of radiative trans¬ 
fer have their origin in the Boltzmann equation for the corre¬ 
sponding particle distribution function J~, in terms of which 

d N = -^T(x,p,t)d s xd 3 p (1) 

is the number of particles within the phase-space volume 
d 3 a;d 3 p, where g is the statistical weight of the species. “Ra¬ 
diation” in the present context is defined as a distribution of 
particles that move with the speed of light c and that are not 
subject to an external force in an inertial frame (p = 0). The 
Boltzmann equation for radiative transfer in a fixed frame 
then becomes (n = p/\p\)\ 

-^-X + n- X x X = B. (2) 

c at 

Here, and in several following cases, we suppress the func¬ 
tional dependencies. The collision integral B = B(x,p,t) 
generally contains explicit integrals in momentum space, 
making Eq. |2| an integro-partial differential equation. In¬ 
stead of working with the distribution function directly, for 
the macroscopic view one prefers using the frame-dependent 
specific (i.e. monochromatic) intensity 

X(x,n,e,t) = (e/hc) 3 cX(x,p,t ), (3) 

wher^] e = \p\c. Bearing in mind that an essential part of 
the collision integral depends on the distributions of fluid 
particles it is often preferable to measure X in the frame 
comoving with the fluid (“comoving frame”, “fluid frame”), 
since in that frame the isotropy of the fluid particle distri¬ 
bution^] induces symmetries in the collision integral that 
make it computationally most feasible. Using arbitrary, but 
fixed, Eulerian spatial coordinates defined in a frame we de¬ 
note as the laboratory frame (“lab-frame”) and momentum 
space coordinates (i.e. e and n) defined in the fluid frame, 
the comoving-frame equation of radiative transfer up to or¬ 
der (D(v/c) (v = |u| is the velocity of the fluid as measured 


2 We will use the terms “energy” and “frequency” interchange¬ 
ably when referring to the corresponding degree of freedom in 
phase space. 

3 We implicitly assume the fluid to be in local thermodynamic 
equilibrium (but not necessarily in equilibrium with neutrinos). 
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in the lab frame) becomes (e.g. Buchler|1979 Kaneko et al. 
1984 Munier fc Weaver|1986 1: 


1 dX 

v ■ n dX 


— 


+ n 

c dt 

c 2 dt 


d 

-r fan 

1 


-r- 

+ - 

de 

V c 

c 


dX v 3 dl 
dx J c dxi 


dn i 


f a ■ n, i a 1 i j k T -7 1 j t—7 i 

—r—n-rl—n n n V,m,- n \7,v 

V c 2 c 2 c c 


r i i k f rri i fc 

j k nn -1 j k vn 


+ X 


an 1. 

r 2 


= c, (4) 


where a = T* fc are the Christoffel symbols associated 
with the spatial coordinates and C = ( e/hc) s cB. Equa¬ 
tion 0 can be derived from Eq. § using Eq. § and the 
0(v/c) versions of the Lorentz transformations for X, t and 

n. 


2.2 Angular moments of the transfer equation 

In order to reduce the dimensionality of the radiative trans¬ 
fer problem and to construct the link to the hydrodynamics 
equations, one utilizes the fact that the specific intensity 
is related to the specific (frequency-integrated) energy den¬ 
sity E ( E ), energy flux density F 1 (F x ) and pressure tensor 
pij (pyj Q f radiation by virtue of its angular moments of 
increasing order, defined by 


{cE, F l , cP XJ , } = / dXlX{l, n, nn J , nWn } , (5) 


and 


{E,F\P x \Q ljk }= I de {E, F\ P lJ , Q lJK } 


/< 


ro'fci 


( 6 ) 


where Q ljk and Q l ^ k are the analog 3rd-moment quantities. 

In the following, we neglect terms including the acceler¬ 
ation a 1 and the second term containing the time derivative 
in Eq. ([4]). These terms are effectively of order 0(v 2 /c 2 ) 
for temporal changes of the velocity and radiation fields oc- 
curing on a fluid timescale l/v, where l is a characteristic 
length scale of changes in the hydrodynamic background 
and v a typical fluid velocity (Mihalas & Mihalas 1984 see, 


however, |Rampp fe Janka 2002 and L owrie et al.||200Tl for 
comments on the second term of Eq. @). Temporal changes 
of these fields on the radiation timescale l /c would enhance 
the importance of the aforementioned terms. In that case, 
however, the preceding validity assumption of the 0(v/c ) 
equation may become questionable to begin with anyway. 
It is worth noting that the results of one test in Sec. |4.3.2| 
suggest that for the conditions in CCSNe the omission of 
the aforementioned terms is justified. 


2.2.1 Moment equations of energy transport 

The system for the first two moments of Eq. © , exclud¬ 
ing the aforementioned terms that are effectively of order 
0(v 2 /c 2 ), is obtained by performing the angular integra¬ 


tions as in Eq. © and it reads 
d t E + VjF j + 

Vj(v j E) + ( V jVk )P jk - (V jVk )dJeP jk ) = C co) , (7a) 
'-.-' v -*-' '-*-' 

(I) (II) (IV) 

d t F i + c 2 X7jP ij + 

VjtfF) + F j VjF - (V J -i»)0«(eQ <J ~ h ) = C (1) ’ i , (7b) 

(I) (HI) (IV) 


where C* 0> = f dSIC and C^’ 1 = f dfl n l C. The labeling 
by Roman numerals denotes the different types of velocity- 
dependent terms: The (I)-terms account for the change of 
the comoving-frame moments owing to advection. The (II)- 
terms account for the change of radiation energy due to com- 
pressional work against the radiation pressure. The (III)- 
terms account for changes of the radiation fluxes due to 
aberration. The energy-coupling (IV)-terms are responsible 
for the change of the spectral shape of the radiation field 
associated with the Doppler shift. Note that in the grey 
formulation of the moment equations (i.e. after integrating 
Eqs. ([Tj) over energy e) the (IV)-terms vanish. For the ex¬ 
plicit form of Eqs. 0 in spherical polar coordinates (addi¬ 
tionally including the aforementioned terms that are omitted 
in our presentation) the reader may consult the Appendix 
of Buras et al. (20061. 

Equations ( 7j) are the radiation evolution equations used 
in our code. Since the source terms in general depend on the 
energy and species of neutrinos, we have to solve a set of 
moment equations for each energy group (after discretizing 
the energy space into a finite set of energy groups/bins, cf. 
Sec. 0 and for each species. Hence, given N t energy bins, 
Nape species and taking into account IVdi m components of the 
flux density F z , we have to process (Adim + 1) x N spe x N e 
equations in total in our multidimensional, multi-group radi¬ 
ation transport scheme. For the following presentation, how¬ 
ever, we will only indicate individual species or the energy 
dependence if it is demanded by the context. 


2.2.2 Moment equations of number transport 

The moments connected to the number transport (number 
density, number flux density etc.) are given by 

{ N, Fj, Pj,Q% k , ...} = e~ 1 { E, F\ P ij , Q ijk ,...}. (8) 

Although we do not directly use them in our code, we list 
the equations describing the neutrino number evolution for 
completeness here. They are structurally similar to Eqs. 0 
except for terms associated with the energy derivatives: 

d t N + VjFj + Vj(v j N) - VjV k dJePj k ) = e-'C^ , (9a) 
d t Fj + c 2 VjPj + VJFF^E 

FjfVjV* - S7jV k [Q% k + dJtQ% k )] = . (9b) 

2.2.3 Transformation into lab-frame 

The transformation of energy-integrated moments from the 
comoving into the lab-frame can be performed by referring 
to their intrinsic tensorial structure which dictates the way 
the Lorentz transformation has to be applied. The energy- 
associated moments F,F l , are components of a 2nd-rank 
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tensor, namely the energy-momentum tensor of radiation, 
while the number-associated Oth and 1st moments combine 
to a 4-vector. This results in the following transformation 
rules correct to order 0{v/c ) for the energy-related mo¬ 
ments: 


Flab = E + 2c 2 VjF° , 

(10a) 

Kb = /• ' + v*E + v,P J , 

(10b) 

>K-P'' • (-• /•' • r’F'). 

(10c) 

and for the number-related moments: 


N\ ab = N + c -2 ViF l N , 

(Ha) 

F i NMb = K+v i N. 

(lib) 


Note that these transformation rules only apply for the grey 
quantities but not for the monochromatic moment^ The 
energy-integrated source terms transform into 

the lab-frame source terms CfabiCjab a similar wa y as 
N and Fn (cf. Eq. (Ill), i.e. as a 4-vector, since they are 
defined to form the right-hand side of a conservation law of 
a 2nd-rank tensor in its original relativistic formulation. 


2.3 Interaction source terms and coupling to 
hydrodynamics 

The interaction source terms are the actual terms that intro¬ 
duce the microphysical properties and the coupling of mat¬ 
ter and radiation into the transport problem. The source 
terms C^°\ C 1 ' 1 ' 1 ’’ 1 for the neutrino moments give rise to cor¬ 
responding hydrodynamic source terms Qm, Q e that ac¬ 
count for the change of fluid momentum and gas internal 
energy, respectively, due to the interaction with neutrinos. 
We restrict ourselves here to the basic Euler equations and 
neglect additional physics, such as viscosity, magnetic fields 
or the co-evolution of a set of nuclear species. The evolu¬ 
tion of the baryonic density p, momentum density pi/, total 
gas-energy density et = e\ + pu 2 /2 (where ei is the internal 
energy density) and electron fraction®] W is then dictated by 


the system 

d t p + X7j(pv j ) = 0, (12a) 

d t (pY e ) + Vj(pY e v j ) = Q N , (12b) 

d t (pv l )+ X7 j (pv t v J + P s ) = Qm, (12c) 

d t e t + Vj (v j (e t + P g )j =Qe+VjQ j m . (12d) 


The gas pressure P g is obtained by invoking an equation of 
state (EOS), which at the same time provides the quantities 
required to compute the opacities (such as the temperature, 
chemical potentials, or particle densities in case of nuclear 
statistical equilibrium). By virtue of the physical meaning 


4 However, corresponding 0(y/c) expressions for the monochro¬ 
matic moments can be formulated in terms of Taylor expansions 
of the moments in energy space (e.g. |Mihalas & Mihalas||1984] 
Hubeny & Burrows|2007 1 . 

5 As usual, this quantity is defined as Y e = (n e ~ — n ( . )/tib 
with the number densities n e ± of electrons and positrons and the 
baryon number density ng. 


of the moments E and F l , the source terms for the hydro- 
dynamic equations can immediately be identified with 


Q E 


Qm 


Q N 



(13a) 

(13b) 

(13c) 


where mn is the atomic mass unit and the sums contain all 
contributions from individual neutrino species. 

At present, the most important types of (electron-) neu¬ 
trino interactions are implemented, namely the capture of 
(anti-) neutrinos and (anti-) electrons by free nucleons and 
nuclei, isoenergetic scattering of (anti-) neutrinos off free 
nucleons and nuclei, and inelastic scattering of (anti-) neu¬ 
trinos off (anti-) electrons. All corresponding source terms 
are adopted as described in Rampp & Jankaj2002l, which 
is mostly based on the compilation given by Bruenn (19851. 

For introducing further concepts in a way that keeps 
the presentation as essential as possible, in the following we 
will for simplicity assume that only isoenergetic scattering 
and absorption/emission reactions take place. In this case 
the source terms in the moment equations, Eqs. (J7|, can be 
written as (e.g. Bruenn 1985]) 


C (0) = c« a (E eq - E) , (14a) 

C (1),i = -c(k a + k s )F 1 , (14b) 


where K a and k s are the combined absorption (corrected for 
final-state Fermion-blocking) and scattering opacities, and 
_E eq is the equilibrium energy density associated with the 
Fermi-Dirac distribution J-fd, 

E^(e,p,,T) = j dn(^) 3 J- FD 

/ € \ 3 1 
~ \hc) exp{(e — /u„)/(fc B T)} + 1 ’ 

which is a function of the fluid temperature T and the chem¬ 
ical potential p u of the corresponding neutrino species. The 
transport opacity Ktra, which is given by K t ra = K a + «s, is 
related to the mean free path A„ between two momentum¬ 
exchanging reactions via Kt ra = A^ 1 . 


2.4 Algebraic closure methods 

The full information contained in the Boltzmann equation 
can be captured equally well by an infinite series of con¬ 
servation equations for the angular moments, in which the 
evolution equation for a moment of rank m contains the mo¬ 
ment of rank m + 1 within the divergence operator. Instead 
of solving the infinite series of moment equations, the series 
can be truncated at the level of the (m + l)-th moment, 
provided the (m+ l)-th moment is available to close the set 
of m equations. However, in order to determine the (m-f-1)- 
th moment that is consistent with the Boltzmann equation, 
it is inevitable to solve the latter in some approximation or 
the other, using computationally expensive methods such as 
discrete ordinate or (long- or short-) characteristics schemes. 
A computationally much cheaper, though more approximate 
option is to assume that an algebraic closure relation holds 
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between the evolved moments and the (m + l)-th moment. 
This is what defines the algebraic closure methods, such as 
FLD and AEF. Essentially, this corresponds to imposing ad¬ 
ditional conditions or symmetries on the local radiation field. 
The consequence is that two (out of seven in the general 
case) independent variables describing the angular depen¬ 
dence of the radiation field disappear from the treatment. 
Evidently, the tradeoff for this computational simplification 
is that an algebraic closure method may strongly vary in 
quality between different physical setups. For example, since 
in an algebraic closure method the consistent evolution of 
higher-order angular moments is ignored, it appears likely 
that the quality of the scheme appreciably depends on the 
geometric complexity of the radiation field, or equivalently 
on the shape and number of the individual radiation sources. 
Moreover, connected to this issue is the circumstance that 
in the optically thin limit of vanishing source terms an al¬ 
gebraic moment scheme is generally not able to accurately 
describe the unperturbed superposition of multiple radia¬ 
tion fronts, simply on account of the closure being a local, 
non-linear function of the evolved quantities. Despite these 
conceptual deficiencies, which have to be individually tested 
for in each case of application, algebraic closure methods 
can in many cases offer an excellent compromise between 
efficiency and accuracy when performing energy-dependent, 
multidimensional RHD simulations. 

Independent of the rank at which the scheme is trun¬ 
cated, any closure prescription should agree with certain 
consistency requirements that directly follow from the def¬ 
inition of the moments or from the Boltzmann equation. 
Using the normalized moments f= F/(cE), where / = |/| 
is the flux-factor, D lJ = /E, the Eddington tensor, and 

qijk — Qijk / {^ C E), the normalized 3rd-moment tensor, it fol¬ 
lows from the definition of the angular moments, Eqs. <©. 
that 


\fi 


i, 

(16a) 

D'3 


i, 

(16b) 

E D jJ 

= 

i, 

(16c) 

\q ijk | 


i, 

(16d) 


= 

Vv-' -E'/"' - /'• 

(16e) 


3 3 3 


must hold at any time. In the “free-streaming limit”, which 
is approached far away from regions of radiation-matter in¬ 
teraction, all of the radiation propagates into a single direc¬ 
tion away from its source and it must hold 

/ = 1 , D lJ = n F n F , , q 3k = n F n J F n F , (17) 

with n F = F l /\F\ denoting the direction of the flux den¬ 
sity. In the opposite limit of very frequent interactions, i.e. 
in the “diffusion limit”, the specific intensity is approxi¬ 
mately isotropic. Ignoring velocity terms, the radiation mo¬ 
ment equations degenerate in this limit to the diffusion equa¬ 
tion 


d t E + V, 


c 


'Tv l ra 



= c (0) 


and the relations 


1 VE 
3Ktra E 



(18) 


(19) 


In the following sections, we will outline the basic concepts of 
FLD and AEF schemes and present several closure prescrip¬ 
tions. In Sec. |4.3.I| we will explicitly compare both methods 
and the presented closures on the basis of the neutrino ra¬ 
diation field produced by a proto-NS. 


2-4-1 Flux-limited diffusion method 


The approach of FLD (e.g. |Wilson et al.|1975||Levermore fe| 
Pomraning 1981) is to truncate the set of moment equations 
at the level of the lst-moment equation and to derive an 
expression for the flux density based on the diffusion limit 


described by Eqs. (181 and (191. Introducing the “Knudsen 


number” R = IJ2I, with 


R = 


X7E 


UJ /‘vtr; 


E 


( 20 ) 


where 


lO = {k-bE + K a E eq )/(KtraE) (21) 

is the “effective albedo”, the flux density Ffld is prescribed 
as 


-Ffld = — A (R) RcE , 


( 22 ) 


in which A (R) is called the “flux-limiter”. The latter is con¬ 
structed such that the flux density correctly preserves the 
constraints of Eqs. |l7|, (1191). To this end the limits 


and 


lim A (R) R = 1, 

R—¥ oo 


lim A (R) = - , 

R —^0 o 


(23) 


(24) 


respectively, have to be fulfilled (note that u> —> 1 in the 
diffusion limit). Three prescriptions are widely used in the 
literature (jWilson et al.|1975 Liebendorfer et al.|2004 Lev 


ermore & Pomraning 19811: 

1 


Av 


JR) = 


3 + R 


ABruenn(R) — 


min ^A wi i son (77), 1+ J X ^! r l _ 


. Av 


JR) 


,r > r„ 
, else, 


ALevermor e(-R) — j-, ^COlll It R^ 


(25a) 

(25b) 

(25c) 


The limiter in Eq. (25b I is only designed for the spherically 


symmetric case in which r is the radius coordinate and r v is 
the radius of a (properly defined) neutrinosphere. The lim¬ 
itation for r > tv is intended to ensure that the neutrino 
flux cannot be higher than if the neutrinos were distributed 
isotropically into a finite cone subtending the sphere of ra¬ 
dius r„. 

The main drawbacks of FLD are: First, the prescription 
of the flux density is in general not consistent with the lst- 
moment equation. As a direct consequence, the full RHD 
system suffers from momentum and therefore energy non¬ 
conservation (Bruenn 1985; Baron et al. 19891 whenever mo¬ 
mentum transfer between matter and radiation takes place. 


In the case of CCSNe, the violation is found (Cernohorsky & 
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A new multidimensional neutrino-hydrodynamics code 7 


van de n HornJT9 90) to be particularly significant in the semi¬ 
transparent region where A (R) undergoes the main part of 
the transition 1/3 —» 0. Even though interim solutions of this 
shortcoming can be formulated, e.g. by introducing an arti¬ 


ficial opacity (Janka 1991 Dgani & Janka 19921 that con¬ 


tains the missing information of the lst-moment equation, 
they introduce further degrees of freedom, hence rendering 
the resulting method rather tuned to special cases. 

Second, in more than one dimension a complication 
arises from the fact that the flux density vector is always 
directed opposite to the gradient of the energy density since 


the pressure is effectively isotropic (cf. Eq. (191): Radiation 
in the free-streaming limit will not keep its original flux di¬ 
rection after closely passing opaque objects. Instead it be¬ 
haves like a gas and fills up space in every direction, unable 
to form persistent shadows. 

The third issue is a purely computational aspect: The 
energy equation evolved in FLD is - at least whenever / 1 

- of parabolic mathematical nature. As such, it comes with 
the property that the operator V ■ F needs to be treated 
time-implicitly in most practical cases. This is because the 
characteristic timescale tfld associated with V ■ F can be¬ 
come extremely short in the optically thin limit K tra —> 0. 
One can roughly estimate the local timescale tfld by think¬ 
ing of the operator V • F as being locally a linear convex com¬ 
bination of the advection operator afcVE and the diffusion 
operator (1 — a)(—Ac/fttra)V 2 E, with some weighting factor 
0 < a < 1. A heuristic dimensional analysis then gives 


a-f—fE + (1 — a)- 

TFLD Ax Ktra Ax 2 


A E 


tfld 


(_+// + (1 ^ q)3A 

\ T’adv 7"diff 


-1 


where Ax is the local grid size and 
T ad v = Ax/C, 

Tdiff — 3KtraAx jc 


(26a) 

(26b) 

(27a) 

(27b) 


are the characteristic timescales of advection and diffusion, 
respectively, ffence, owing to the fact that Taur —> 0 in opti¬ 
cally thin regions, tfld can drop significantly below r a dv 


2.f.2 Algebraic Eddington factor method 

In the AEF method the truncation of the moment equations 
takes place at the level of the 2nd-moment equation, i.e. the 
Eddington tensor is expressed as a function of evolved quan¬ 
tities. Besides resolving by construction the first two draw¬ 
backs of the FLD scheme mentioned in Sec. |2.4.1[ another 
important computational difference to the FLD scheme is 
that the equations solved in the AEF scheme are intrinsi¬ 
cally hyperbolic, which allows to use explicit time integra¬ 
tion methods (at least for all but the source terms, see Sec.[3| 
with time steps not lower than the minimum of r a dv over the 
computational domain. 

Historically, algebraic closures have been constructed 
and analyzed most often in the context of one-dimensional 
systems and up to the present day quite a number of 
one-dimensional closures have been proposed in the liter¬ 
ature. The algebraic Eddington factor x = P/E is typi¬ 
cally expressed as x = x( e i /)> where e = (47r) _1 f dflF = 
(he) 3 /(4ne) 3 E. Note that in constrast to the FLD scheme, 


the closure only depends on radiation moments and not on 
the opacities. In the course of this paper we consider a vari¬ 
ety of closures, of which the Eddington factors are given by 


XMinerbo (/) = \ (6 f ~ 2 f + 6 f) 


XM 1 (/) = 


3 15 

3 + 4 f 
5+ 2^4-3 f 2 


X Ja n ka (/) = Ul + i / 1 ' 31 + ^/ 3 ' 56 


(28a) 

(28b) 

(28c) 


XMax.Ent . (e, /) = |(1 - e)(l - 2e)cr (jzTt) + \ ’ ( 28d ) 


where a(x) = x 2 (3 — x + 3x 2 )/5 in Eq. (28d|. The statistical 
closure XMinerbo(/) by Minerbo) ( j!9781 assumes a fermionic 
particle distribution with low density (e —¥ 0) to be in the 
state of maximum entropy. The generalization of this clo¬ 
sure is XM a x.Ent.(e, /), derived by Cernohorsky & Bludman 


(1994), which additionally to XMinerbo takes into account 
the effects of fermion blocking for e > 0. In Eqs. (28aI, 


(28d| we employed the polynomial approximation derived 


by Cernohorsky & Bludman (1994) to circumvent the nu¬ 


merical inversion of the Langevin function occurring in the 
original formulations of both closures. The Mi closure can 
either be derived from the assumption of the radiation field 


being isotropic in some unspecified frame of reference (Lev- 


1984 I or from maximizing a photon entropy func¬ 
tional (Dubroca fe Feugeasfl9991. Note that in both afore¬ 


mentioned derivations the M i closure actually relates the 
energy-integrated moments. In our present treatment, in 


contrast, we apply all closure relations in Eqs. (281 indi¬ 


vidually for each neutrino energy e. Finally, the closure by 


Janka (1991) was obtained by fitting the neutrino radia¬ 


tion profile around a proto-NS calculated with Monte-Carlo 
methods. For further detailed discussions about specific one¬ 
dimensional closures and their properties, see, e.g., |Smi 
et al. (20001 and |Pons et al. (20001. 


3 


To extend the one-dimensional Eddington factor x to 
the multidimensional Eddington tensor D lJ and the 3rd- 
moment tensor q 13 , we make use of the underlying assump¬ 
tion that these tensors are only a function of the scalar E 
and the vector f] From symmetry arguments it follows that 
D lJ must be symmetric with respect to rotation around the 
flux direction, rip = F/\F\, which is only fulfilled if D 13 is 
a linear combination of the two tensors <and n F n F (e.g. 


P ennisi fe Trovato||l987 ). After using the trace condition of 
Eq. ( |16c| ) the two remaining coefficients can be expressed as 
functions of a single parameter, X: which is the multidimen¬ 
sional generalization of the Eddington factor and is defined 


f dfl (n ■ n F ) 2 T 
X = fdflF 


(29) 


The Eddington tensor then reads (e.g. Levermore 19841: 


D ij = 1 — X 5 ij + — 1 n/ F n j F 


(30) 


6 Mathematically speaking, we assume D^ and q^ k to be 
isotropic tensor functions (e.g. |Pennisi &; Trovato||1987| ) of E 
and F, which in particular implies that these functions may not 
depend on derivatives of E or F. 
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8 Just, Obergaulinger, & Janka 


The energy-dependent comoving-frame moment equa¬ 
tions, Eqs. |7|, also contain the 3rd moments, q 1 ^. In anal¬ 
ogy to the derivation of D above, the condition that this 
tensor only depends on E and F must result in be¬ 
ing invariant under rotation around n F , which is only ful¬ 
filled if q llk is a linear co mbination of re i F <5 J ' fc , nj,5 lk and UpS™ 


(e-g- 


Pennisi 


as well as 

coefficients can be eliminated using the trace conditions of 


1992). The corresponding 


Eq. (16eI in favor of a single parameter, q , defined as 


f dfi (n ■ uf ) 3 T 
q= f dH T ' 
yielding for the 3rd-moment tensor: 


(31) 


ijk _ 


f-q 


(n F 5 jk + n j F 8 ik + n k F S ij ) + 


5q-3f 


■ n F n J F n F . 

(32) 

The 3rd-moment factor q explicitly depends on the 
distribution function. Therefore, only closures that dictate 
an explicit functional form of the distribution function are 
suited for the computation of the 3rd moment, unless ad¬ 
ditional assumptions are made in the construction of the 
closure. For the Minerbo closure, the factor q can be calcu¬ 
lated in a straightforward manner in analogy to the deriva¬ 
tion of x I Minerbop978 1 and reads (again using the poly¬ 
nomial approximation of the Langevin function as given in 
Cernohorsky fe BIudman[l9941 


9Minerbo(/) = ^ (45+10/-12/ 2 -12/ 3 +38/ 4 -12/ 5 +18/ 6 ). 

(33) 

We outline the derivation of Eq. (331 in Appendix [A] Re¬ 
cently, Vaytet et al. (2011) calculated the 3rd-moment factor 


q also for the M\ closure. 


3 NUMERICAL METHOD 

3.1 Motivation of the integration scheme 

Before presenting the details of the numerical implementa¬ 
tion of the AEF scheme outlined in the previous section, we 
briefly summarize and motivate the general framework of the 
discretization scheme. Owing to the fact that the advection- 
type operators on the left-hand side of the two-moment sys¬ 
tem, Eqs. 0 . are of hyperbolic mathematical nature, we 
can employ Godunov-type finite-volume methods commonly 
used in numerical hydrodynamics to discretize these oper¬ 
ators. However, in regions of strong coupling with matter 
the source terms become stiff and the moment equations 
approach the parabolic diffusion limit. Hence, the time inte¬ 
gration is performed in a mixed explicit-implicit manner, in 
which all terms on the left-hand side of Eqs. 0 are treated 
explicitly in time while the source terms on the right-hand 
side of Eqs. 0 are handled implicitly (at least whenever 
being in the stiff regime). In that way the overall time step 
used for integration of the whole scheme is constrained by 
the Courant condition to be on the order of the advection 
timescale r a d v = A x/c, i.e. the light-crossing time of grid 
cells with width Ax. The alternative would be to integrate 
the full two-moment system implicitly in time. In that case 
the computational cost per time step would be significantly 
higher (particularly in the multidimensional case) but on 
the other hand this would allow to employ a larger time 


step which is closer to the fluid-dynamical time step Ax/v. 
We opted for the former version of integration, mainly for 
the following reasons: 

(1) Since the velocities in CCSNe and NS mergers are 
typically rather high, v ~ 0(0.lc), the characteristic hydro¬ 
dynamics time step and therefore the implicit radiation time 
step turn out to be only a factor of a few greater than the 
explicit radiation time step. (2) Since all operators contain¬ 
ing spatial derivatives are treated explicitly, the common 
parallelization methods can be applied with very high effi¬ 
ciency. This is particularly advantageous in the multidimen¬ 
sional case in which a fully implicit scheme would become 
increasingly expensive (as its computational cost typically 
increases faster than linear with the number of grid zones). 
(3) Light fronts and discontinuities in the radiation field 
can be sharply resolved, which tend to be smeared out in an 
implicit method, unless a time step comparable to the ex¬ 
plicit one is used. (4) The overall numerical implementation 
is less involved than for an implicit scheme because inver¬ 
sions of large matrices that couple spatial grid points are 
avoided. (5) While high-order spatial reconstruction meth¬ 
ods can be implemented in a straightforward manner in ex¬ 
plicit schemes, they are usually too prohibitive to be used 
in implicit schemes. 

For other implementations of a similar explicit-implicit 


integration method see,e.g. Sadowski et al. (2013), Skinner 


& Ostriker (20131, O’Connor (20141. For implementations 


of fully implicit schemes see, e.g., Audit et al. (20021, Hayes 


& Norman (2003). 


3.2 Basic discretization scheme 

The spatial and energy-space discretization scheme for all 
quantities is based on the finite-volume approach. The set of 
spatial coordinates can be varied between Cartesian (x, y, z), 
cylindrical and spherical polar (r, 9, <j >) coordinates. How¬ 
ever, for outlining the discretization scheme in this section 
we will restrict ourselves to Cartesian coordinates. The ex¬ 
tension to curvilinear coordinates is realized by adding the 
appropriate geometric source terms to the presented dis¬ 
cretized derivatives. The geometric source terms are purely 
algebraic functions of the evolved quantities and the grid co¬ 
ordinates and are discretized simply by replacing the argu¬ 
ments of these functions with the corresponding discretized 
quantities. 

The spatial grid is composed of volume cells (i, j, k) that 
are obtained after discretizing a given domain in each coor¬ 
dinate direction {x,y,z} into {N x , N y , N z } zones, of which 
each is defined by the cell boundaries * i± i, y i± i, z k± i, with 
{i, j, k} = {1... N x , 1 • • • N y , 1... N z }. The cell-center coor¬ 
dinates are computed as x\ = l/2(x i _ i + x- t+ i ) and analog 
for the other directions. The volume of cell (i, j, k) is denoted 
by AUj.k, and the area of the surface (i + ^,j, k), located be¬ 
tween cells (i, j, k) and (i + 1, j, k), is denoted by AA i+ 1 . k . 

For the grid in energy space, given by N e energy bins, 
we use i to denote the boundaries of the £’th bin, with 
£ = 1... Ne. Furthermore, ej = l/2(e^ + i +e 4 _ i) and Ae^ = 
e, i — e, i define the center and width of the f’th bin, 

s“T 2^2 

respectively. 

We define the radiation fields as well as the fluid quan¬ 
tities on the same spatial grid. The discrete representa- 
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A new multidimensional neutrino-hydrodynamics code 9 


tions U € {p, {pY e ), (pv 1 ), e t } of the hydrodynamic quan¬ 
tities U £ {p, pY e , pv l ,e t} are interpreted as cell-volume av¬ 
erages in space, i.e. as 


t/ij,k = 


AM. 


,j,k J AVj ; 


dV U , 


(34) 


while the discrete representations X £ {E. F\ P IJ , 
of the radiation moments A' £ {E, F l , P l \ are inter¬ 

preted as averages in space and integrals in energy space in 
the following way 


*u,k,e - 


l 



deX . 


(35) 


The spatial and energy discretization of the moment equa¬ 
tions, Eqs 0’ is realized by applying the spatial averaging 
and the energy integration as in Eq. ( |35| > to the moment 
equations. 

The discretization of the spatial derivatives demands it 
to reconstruct quantities from cell-volume averages to call- 
face averages, which for example on the cell-face (i + j, k) 
are given by 


X, 


l — 


+Ti,k,e 


A A; 


/ dA f 

J AA. , 1 . , J At 


deX. 


(36) 


+ i j,k ma , 


The reconstruction algorithms that we use are adopted from 
the hydrodynamics part of the code and can be switched 
between piecewise-constant, piecewise-linear and high-order 


“monotonicity preserving” (MP) schemes (Suresh & Huynh 


1997) or “weighted essentially non-oscillatory” (WENO) 
schemes ( Liu et ah]|1994| |. In what follows, we symbolically 
use X L and A' K to denote the reconstructed values of a 
quantity A' on the left- and right-hand side of an interface, 
respectively. 


3.3 Summary of the integration algorithm 

In order to formulate the algorithm to integrate the evolu¬ 
tion equations for the radiation moments, X = (E , F) (cf. 
Eqs. 0), and for the fluid quantities, U = (p, pY e , pv, e t) 
(cf. Eqs. ( |12| )), we decompose these equations in the follow¬ 
ing way: 


d t X + (<5 t X) hyp + (S t X) vel — (^t X). 

dtU + (<5tU)h y d = («tU) B 


(37a) 

(37b) 


In Eq. (37a), (5 t X) hy p = (VjF- 7 ', c 2 VjP IJ ) represents the 
velocity-independent hyperbolic advection terms, while all 
velocity-dependent terms are subsumed in (d t X) ve i, and 
the interaction source terms are represented by (5tX) src . 
In Eq. (37b l, (<5tU) 3r c are the radiative source terms given 
by Eqs. (131, while all remaining terms representing non- 


radiative physics are contained in (<5iU)h yd - We first sum¬ 
marize the overall integration scheme of Eqs. (371. Subse¬ 


quently, in Secs. [3~4}{3.7| we will explicitly describe how the 
individual terms are computed. 

The following steps are performed to evolve the RHD 
equations, Eqs. (371, from time t n to f n+1 = t" + At. Note 


that we use the shorthand notation A" = A(t n ) to refer to a 
quantity A at some time step f n : 


1) Compute the global integration time step At used for 


© 2015 RAS, MNRAS 000, |T]{30] 


both hydrodynamics and radiation transport as (i £ 


A t rad = min 


( A *)i,j,k 


i,j,M [ |*4kl + max|(A±)L k | 

(A*)E k 


At hyd = min . 

MM I max |(Aflui d )f Jik 

At = CFL ■ min (At ra d, At hy d} 


(38a) 

(38b) 

(38c) 


where A± and Afl u id are the characteristic velocities of 


the radiation system (cf. Eqs. (46 47)) and of the fluid, 
respectively, CFL is the chosen Courant-Friedrichs-Lewy 
number and (A*)T k is the length of cell (i, j, k) in coor¬ 
dinate direction i. 

2) Construct the advection operator (<5 t X)" dv = (<S t X)£ yp + 
(<5tX)" el as a function of the radiation moments X" and 
fluid quantities U n using the discretization rules de¬ 
scribed in Secs. 13.41 and 13.51 

3) Compute the hydrodynamics evolution operator 
(<5tU)h yd as a function of U n . Eventually, add other 
properly discretized terms to (<5tU) dyd corresponding to 
additional, non-radiative physics (such as gravitation or 
magnetic fields). 

4) Perform an intermediate update to the time f n+ 5 = t" + 
At/2 by solving 


X n+ s =X n + y 


-(StX)l dv + (S t X) 


n,n+2 

src 


U n +2 = u n + 


At 


-(fctOEyd + (<5 t U)s(c +5 


(39a) 

(39b) 


for X n+ 2 i U n+ 2 ; where the two comma-separated su¬ 
perscripts indicate that the source terms can generally 
depend on the hydro- and radiation variables at both 
the old and the new time step, accounting for the fact 
that the integration can be performed explicitly and/or 
implicitly as the circumstances require (cf. Sec. 3.6 1 . 


5) In analogy to steps 2) and 3) use X n+ 2 and U 1 
compute (<5 f X)^ v 2 and (<5tU)" lyd , respectively. 

6) The quantities at t n+1 are finally obtained by solving 


to 


U 


n + l 


X" + At 


= U n + At 


-(«.*)2? + (**)S i ^ 1 , (40a) 


~(StV)ZI + (&U) 


i+»,n + l 


(40b) 


for X n+1 , U n+1 in an analog manner as in step 4)■ 

The above update scheme is formally unsplilQ and is 2nd- 
order accurate in time with respect to all explicit (advec¬ 
tion or source) terms and lst-order accurate with respect to 
implicit source terms. Although realizations of higher-order 


implicit-explicit (IMEX) schemes exist (e.g. Ascher et al. 
1997 Pareschi & Russo 2005 McKinney et al. 20141 we 


7 By unsplit , we refer to the property that the advection and 
source terms are integrated within a single step, in contrast to 
which in a split scheme the quantities would be updated first 
for one set of terms (including the recomputation of primitive 
variables such as temperatures, opacities and Eddington factors) 
before calculating the other set of terms. 




















































10 Just, Obergaulinger, & Janka 


found that the method described above is sufficiently robust 
and accurate in all applications we considered so far. 


3.4 Velocity-independent hyperbolic part 

Our basic treatment of the velocity-independent, hyperbolic 
part of the two-moment system, Eqs. 17], follows along the 


lines of Pons et al. (20001 and Audit et al.|( 20021. The notion 
is to exploit a Godunov method ( Godunov|19 59) as the ba¬ 
sis for a high-resolution shock capturing scheme that solves 
the local Riemann problems between discontinuous states at 
the interfaces of grid cells. We start the presentation of its 
working method by considering the one-dimensional system 


d t 


+ d x 


F 

c\E 


= 0, 


(41) 


where the algebraic closure X = x( e > /) is a function of e 
and /. This system is hyperbolic if the Jacobian matrix J 
of the vector of fluxes (F,c 2 xE) T , 


J = 


0 1 

c 2 (x + e§?-/§f) c§* 


df 


(42) 


has real eigenvalues A± , given by 


.id cdx.c dx 2 ,9\ f 9x, 

X± =2df ^df +4(x + e dk- f df } 


(43) 


All of the closures listed in Eqs. ( |28| ) fulfill the condition 
of hyperbolicity and lead to the following properties: In the 
free-streaming limit, / —> 1, we have 


i \ 1 D . \ 1 D ,d X 

X - 1 ! — + C ! — ( ~qJ — 1)C, 

while in the diffusion limit, / —> 0, one obtains 

1 ,id , 1 

X = - , A + = ±—=c. 

3 ^3 


(44) 


(45) 


That is, the limiting cases for the Eddington factor and 
the wave speeds are consistent with what is dictated by the 
Boltzmann equation. 


In the multidimensional generalization of Eq. (411 the 


matrix eigenvalues contain an additional dependence on the 
direction cosine y = cos cxf, where of is the angle between 
the direction of the radiation flux vector F and the coordi¬ 
nate direction with respect to which the derivative is taken. 
The wave speeds are now obtained as roots of a cubic polyno¬ 
mial leading, at least in terms of a general closure, to rather 
large expression^ For practical purposes we do not take 
into account the exact angular dependence of the eigenval¬ 
ue^] A™ act (fi) but we instead approximate the latter using 
the following lst-order expansion in y. 


a±( m ) = xt 


fc (0) + \fi\ (A± xact 


(1) - A± xact (0)) , (46) 


8 See, however, |Skinner &; Qstriker||2013| who found for the par- 
ticular case of the M\ closure comparably compact expressions for 
the wave speeds as functions of p, and /. 

9 Note that also a third eigenvalue Ao appears in the multidimen¬ 
sional case which fulfills A_ < Ao < A+. However, this eigenvalue 
is not relevant for our present purpose. 


where 


\ exact / -j \ \ ID 

A ± UJ — A ± 


\ exact 
A ± 


(°) = 


(47a) 

2(l-X-e|) + i0(l + 2/ 2 -3 X ). 

(47b) 


In Appendix[B]we provide the components of the Jacobian in 
terms of a general closure x(e, /) and show plots of the exact 
and linearized wave speeds for some specific closures. These 
plots indicate that the linearized wave speeds reproduce the 
exact wave speeds sufficiently well for the former to be used 
instead of the latter as estimates for signal speeds of an ap¬ 
proximate Riemann solver (see below). It is worth to note 
that the qualitative behavior of the angular dependence of 
the wave speeds, A±, is physically consistent with the under¬ 
lying Boltzmann equation: In the diffusion regime, / 1, 

the wave speeds become nearly isotropic with |A±| —> c/y/3, 
while in the free-streaming regime, f —> 1, the wave speeds 
become forward-peaked with A± — > c in the direction of F 
and A± —» 0 orthogonally to F. 

In a fashion that is commonly employed in numerical 
hydrodynamics, we use the above velocities as signal speeds 
for an approximate Riemann solver in order to compute 
the numerical fluxes through each cell interface. We use 
the two-wave solver by Harten, Lax and van Leer (HLL, 
Harten et al. 1983], which approximates the final numeri¬ 
cal interface fluxes as functions of the left-/right-hand side 
fluxes F l / r (with F £ {F l , c 2 P lJ }) and states U L/,R (with 
U £ {E,F*}) as 


\JrinC17C .HL 
pHLL _ A + r — 


XHLL _ X HLL 


+ 


A HLL A HLL (yR _ yL) 

a hll _ x HLL 


(48) 


with the signal speeds A RLL = max(0, A+, A r ) and A RLL = 
min(0, A r , A r ). All quantities labeled by L/R in this flux 
formula are computed using the cell-interface reconstructed 
moments E L / R , J? n ’ L / R . Applying this solver, the final 
spatially-discretized version of the operator (JfA')hyp of 
Eq. (37a) reads (using X £ {E,F}) 


0*.i,m) 


AA.i . k F RL i L ., , - A A. 1 . ,F HL i L . , , 

2 ’-I’* l +2 , J’k , £ 2 ’J’k l— o’J>k’£ 


hyp 


AK 


U,k 


+ V + u z” 


(49) 


where we symbolically indicated the contributions from the 
y- and 2 -directions, which are obtained in an analog manner. 

Yet, there is a caveat we have to face when approaching 
the parabolic diffusion limit (cf. Eq. ©) with the scheme 
described above, since the latter is originally designed only 
for hyperbolic systems. In contrast to the hyperbolic sys¬ 
tem, the parabolic diffusion equation is not associated with 
characteristic waves propagating information between cells 
with finite speeds. Hence, the ansatz of using a Riemann 
solver that tracks characteristics via upwinding and cap¬ 
tures shocks by adding diffusivity is no longer justified in 
the parabolic diffusion regime. Instead, the fluxes in the dif¬ 
fusion regime should be of central type (i.e. symmetric with 
respect to the cell interface) and they should lead to as lit¬ 
tle as possible numerical diffusivity in order not to spoil the 
effects of the physical diffusivity. To handle this issue, we 
employ a simple switch between the two types of fluxes ac- 
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cording to: 


pHLL,* 

'+h 




(50) 


where the index “i” denotes a representative grid index for 
any coordinate direction and the “stiffness parameter” 


V = fttraA* = 


Ax 

X7’ 


(51) 


is a measure of the degree of neutrino-matter coupling rel¬ 
ative to numerically resolved scales of length and time: For 
V )> 1 neutrino interactions proceed on spatial and temporal 
scales smaller than the grid scale Ax and shorter than the 
numerical time step Ax/c, respectively. Hence, for V exceed¬ 
ing unity the source terms become stiff and thereby cause 
the two-moment system to undergo the transition from a hy¬ 
perbolic to a parabolic system. Our experience from several 
tests (cf. Sec. [4] particularly Sec. |4T2| ) has shown that the 
seemingly discontinuous jump between both flux formula¬ 
tions in Eq. ( |50| ) has no significant influence on the solution. 
This is because at the point of transition, V = 1, one often 
has the situation that (1) the flux factors / are sufficiently 
small to lead to nearly equal contributions of F L and F R in 
F hll , and (2) the relative importance of the diffusive part 
of F iill (i.e. the second term in Eq. (481) is still negligible, 
particularly when using high-order spatial reconstruction. 


3.5 Velocity-dependent terms 

In the following we present the recipes used to discretize 
the velocity-dependent terms (5t-Y) ve i. In order to discretize 
the terms containing velocity derivatives, we reconstruct the 
velocities to obtain {ffi L / R located at each cell interface by 
using the same reconstruction algorithm as for the radiation 
moments. 


3.5.1 Fluid-advection terms 


As fluid-advection terms we denote the (I)-terms in Eqs. ([7|. 
For their discretization we also employ an HLL-type Rie- 
mann solver for each coordinate direction analog to Eqs. ( |48[ ) 
and (491. Specifically, we first compute fluxes F HLL,adv de¬ 


fined by the right-hand side of Eq. (481, but with the nu¬ 
merical interface fluxes F L/,R = fi L / R X lj/R anc j the signal ve¬ 
locities A RLL = max(0, v L , v R ) and A RLL = min(0, fi L , D R ), 
where the v L ' R are the reconstructed velocity components 
normal to the interface at which the corresponding numer¬ 
ical flux is computed, and A' L / R € {E R ^ R , / r ‘> L / R j are the 
reconstructed moments defined at the same interface. The 
final fluid-advection terms are then discretized exemplarily 
in ^-direction as 


M,j,k,f) 

\ / £ 


A a i-HLL,adv a a i-HLL,adv 

AA., 1 . k K i ’ — AA_1 j k h 1 . 

i+ 2>J’ R l+5J,k,i 1 2’J' k l-jj,k,£ 


AHj.k 


and analogously in the other coordinate directions. 


(52) 
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3.5.2 Velocity derivatives 


The remaining velocity derivatives are discretized exemplar¬ 
ily in ^-direction as 


d x v 



f ,j,k 


— A A 


AK 


J,k 



(53) 


while to obtain unique interface velocities u* i ., we arith- 

1 ' 2 ’J ,K 

metically average the reconstructed velocities: 


;,j,k 


^(<4j, k +<; R Uk ) ■ 


(54) 


The discretization of the remaining components of djV x is 
given by analog expressions. 


3.5.3 Doppler shift terms 


In our multi-group treatment of comoving-frame radiation 
transport, we allow radiation energy to be redistributed be¬ 
tween energy groups via the (IV)-terms in Eqs. 0 describ¬ 
ing Doppler shift. From the computational point of view an 
important property of the Doppler shift terms is that they 
have a different functional structure in the energy-based mo¬ 
ment equations, Eqs. |t]|, than in the number-based moment 
equations, Eqs. albeit being physically equivalent. As a 
consequence, a naive discretization of the Doppler terms in 
the energy equation will generally lead to non-conservation 
of neutrino number and therefore of lepton number. Al¬ 
though the non-conservation could be avoided by solving the 
number-based moment equations in addition to the energy- 


based versions (e.g. Rampp & Janka 20021, this would at 


least double the computational expense. We therefore im¬ 
plemented the number-conservative method developed by 


Muller et al. (2010). For a detailed description of this scheme 


we refer the reader to the original paper; in the following we 
only briefly summarize the key features. 

Suppressing spatial grid and tensor indices, we write the 
combined Doppler shift terms of the Oth-moment equation 
for the £’th energy bin as 


(SA) 


Doppler 


= w(P f +F 5 _1 -F ?+ i), 


(55) 


where w subsumes the discretized velocity derivatives, 
denotes a discretized component of the 2nd-moment tensor 
obtained by applying the closure relation to the discretized 
moments and F^, and F ?± i are the discretized versions 
of the effective fluxes eP located at the energy-bin inter¬ 
faces e^ ± i. With the constraint that the energy-integrated 
number density shall be conserved, i.e. that 


E e « _1 (^) 


Doppler 


= o, 


(56) 


the interface fluxes F^ ± i can be written a^j 


F ?+1 = F r +F r 


«+i i 


(57) 


10 Note that these interface fluxes are not uniquely determined 
but chosen as a specific set fulfilling the imposed constraint of 
total number conservation. 
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with 


Ft = 


1 - 


e 5 £ «+i 


F? = 


£ e<7-i 


_i _ 1 -F’c( 1 -7e)- 


(58) 

The weighting factor 75 € [0,1] can be chosen manually 
and we fix it similarly as in |MiilIer et al. (2010). At the 
lower boundary in energy space we usually have the min¬ 
imum energy = 0 and therefore F 1 _i = 0. For the 

upper boundary at £jv £ + i we either use an exponentially ex¬ 
trapolated energy distribution for the high-energy tail or the 
condition that the numerical flux vanishes. For the energy 
derivative of the 3rd moments occurring in the lst-moment 


equation, Eq. (7b I, we use an analog prescription as the one 


given above for the 2nd moments. Specifically, the Doppler 
shift terms for one component F% of the lst-moment vector 
(again suppressing tensor indices) is discretized as 

(StPA =w f (-^)de 

\ /Doppler V de ) 

= w(F e _i -F e+ i), (59) 

where w again represents velocity derivatives and F £ _ 1 are 


calculated via Eqs. (571 and (58 1 exactly as explained above 


but with Pj replaced by Qj. 


3.6 Interaction source terms 

As already mentioned before we intend to use a time step 
At that is close to the radiative advection timescale T a dv = 
Ax/c to integrate the full system of moment equations. 
However, the numerical integration of the interaction source 
terms (J t X) arc , Eq. (37aI, deserves special care because the 


characteristic neutrino-interaction timescale, 
Tint = (cKtra) — A1//C , 


(60) 


can become shorter than r a dv by up to many orders of mag¬ 
nitude. In this case, i.e. for stiffness parameters V > 1 (cf. 
Eq. @), the moment equations are stiff and a fully ex¬ 
plicit time integration would lead to numerical instability. 
Hence, the source terms make an implicit treatment indis¬ 
pensable. However, it is important to note that the charac¬ 
teristic timescales of the advection-type terms (St X)h yP and 
(<5tA') ve i (cf. Eq. (37a)) are longer than r a d v also for V > 1: 
In the diffusion limit the characteristic timescale of (<5 t X)h yP 
is the diffusion timescale Tdiff = 3/t t r a Aa; 2 /c = 3 V ■ r a dv > 
T a dv, while the characteristic timescale of (<5tX) ve i is in¬ 
dependent of the stiffness parameter and roughly approx¬ 
imated Inp^l ~ Ax/v. Hence, we can safely apply a scheme 
in which the operators (<5tX)h yp and (JtX) ve i are treated 
explicitly, while the local contributions from the interaction 
source terms are integrated implicitly in time. Compared to 
an implicit treatment of the full system, this mixed-type in¬ 
tegration greatly reduces the size of the matrix that needs 


11 This estimate disregards fluxes in energy space mediated by 
the Doppler shift terms (see Sec. |3.5.3| l, which in principle can 
cause these terms to change on timescales shorter than Ax/v. In 
practice, however, these timescales are usually longer than r a( j v 
such that an explicit integration of the Doppler terms with time 
step ~ r a ,] v turns out to be unproblematic. 


to be inverted since all spatial derivatives that couple neigh¬ 
bouring cells on the spatial grid are handled by the explicit 
part. 

Depending on the included types of interactions, the 
source terms (JtX) arc and (<5tU) S rc can in general each de¬ 
pend on all evolved radiation and hydrodynamic quantities. 
However, (standard) neutrino interactions do not change the 
baryonic mass density p, and they only have a marginal im¬ 
pact on the fluid momenta pv , at least in typical situations 
where neutrino transport is relevant; hence these quantities 
may be treated explicitly in time. Still, if we were to in¬ 
tegrate all but the aforementioned quantities implicitly - 
which means expressing the source terms (<5tX) src , (<5 t U) arc 
as functions of these variables defined at the new time step in 


Eqs. (391 and (401 - we would generally need to solve a non¬ 
linear system of equations of rank (iVdim + 1) x N c x N spe + 2 
(recalling that A/iim and N spe are the number of evolved lst- 
moment components and neutrino species, respecively, and 
N e is the number of energy bins), in which all radiation 
moments E, F as well as the gas-energy density e\ and the 
electron-number density n e = pY e /m b are handled implic¬ 
itly. However, an implicit treatment of all of these quanti¬ 
ties is not always necessary. That is, under certain condi¬ 
tions a significant reduction of computational expense can 
be achieved by treating a subset of variables explicitly, i.e. 
by using the quantities defined at the old time step in the up¬ 


date formulae, Eqs. (391 and (401. Below we list the different 


modes of the source-term treatment, which we implemented 
in order to avoid a fully implicit integration whenever this 
appears justified: 

a) All radiation moments E, F plus the hydrodynamic vari¬ 
ables Si,n e appearing in the source terms (StX) src and 
(5tU) arc are defined at f n+1 , i.e. the systems of Eqs. (|39|) 


and (401 are solved fully implicitly. Consistently, also 


most of the primitive variables (temperature, opacities 
etc.) are handled implicitly. Only the normalized 2nd- 
and 3rd-radiation moments x and q (cf. Sec. |2.4.2 1 , re¬ 
spectively, and the Legendre-coefficient matrices for re¬ 
actions coupling multiple energy bins (see, e.g., Rampp 
& Janka 20021 are taken from the old time step. 

b) Like a), but the gas-energy density Si and electron- 
number density h e are taken from the old time step. 
This reduces not only the dimensionality of the coupled 
non-linear system by 2 , but also alleviates the compu¬ 
tational expense by the demands of re-computing the 
temperatures and opacities within each iteration step of 
the root finding procedure. 

c) Like b), but all energy-coupling interactions (e.g. 
neutrino-electron scattering) are treated explicitly in 
time. The remaining source terms corresponding to 
emission/absorption and isoenergetic scattering can 
then be written as in Eqs. (141, which results in the 


source terms to completely decouple from each other, 
allowing for a straightforward implicit integration with¬ 
out any matrix inversion. 

Since the criteria for selecting a certain integration mode 
are highly problem dependent, we do not specify them here; 
see, e.g., Sec. |4.3.2| for a particular choice in the CCSN 
context. In any case where a coupled non-linear system of 
equations has to be solved, we make use of the routine 
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nag_nlin_sys_sol from the NAG librarwhich employs the 
“Broyden method” for root finding. 

Based on our experience (see, e.g., last paragraph in 
Sec. |4.3.2| ) and in agreement with O’Co nnor] (] 2014), the 
use of integration mode a) (i.e. an implicit treatment of 
Si and h e ) is barely necessary in practice, even in regions 
of very high stiffness parameters, P > 1. This is because 
in most situations where V 3> 1 neutrinos are trapped in 
the gas and very close to weak equilibrium, which means 
that the net (absorption minus emission) source terms for 
Si and he are effectively small. Consequently, under these 
conditions e\ and n e essentially change only on grounds of 
fluid-dynamical effects on fluid timescales, the latter being 
longer than the radiative timescale r a dv ~ A x/c used for 
time integration. Conversely, an explicit treatment of ei and 
he would not be feasible if we would integrate the radiation 
moment equations fully implicitly using a time step longer 
than ~ Ar/t) (as it is done, for instance, in most existing 
neutrino-hydrodynamics codes using FLD or a Boltzmann- 
solver). Additionally, it is worth to note that although an 
explicit treatment of ei and h e may slightly reduce the ac¬ 
curacy it does not significantly harm the numerical stabil¬ 
ity of the overall scheme, because the effective source terms 
for ei and n e are computed via Eqs. (131 using the time- 


discretized source terms for E v ^, which themselves result 
from an implicit (and hence stable) integration in E v ,£. 

Once after appropriately time-discretized expressions 
for the source terms of the 1st radiation moments are found, 
the changes of the fluid-momentum and kinetic-energy den¬ 
sities due to momentum transfer with neutrinos obtained 
using Eq. |T3b| as: 



(61a) 


(61b) 


Finally, since this is a non-trivial and important as¬ 
pect, we now demonstrate that the time-integration algo¬ 
rithm presented in Sec. |3.3| in combination with the stiffness- 
parameter dependent flux formulation, Eq. ( |50| ), allows the 
radiation flux F in diffusive regions (in which V 1) to re¬ 
lax to the corresponding diffusive flux Fn = — cVf?/(3Ktra) 
in a numerically stable and non-oscillatory manner. To this 
end, we may neglect velocity-dependent terms - which are 
reduced by a factor of 0{{v / c){l /V)\ compared to the domi¬ 
nant terms in our comoving-frame formulation and therefore 
subdominant in the diffusion regime - and we only consider 
the (usually dominant) emission-/absorption and isoener- 
getic scattering reactions, for which the lst-moment source 
terms can be written as ( 5t.F) S rc = —cKtra-F. Since in the 
diffusion regime we expect the pressure tensor to be al¬ 
most isotropic, ~ S 13 E/3, the hyperbolic operator for 
the flux density then reads exemplarily in ^-direction (cf. 
Eqs. j49|, pf): 


(<5tK 


x\ n 
/hyp 


=* -C«tra(FS)r 


E h i + s i+1 

1 ' 9 1 I O 


- A A; 


1 + B -_ 


(62) 


www.nag.co.uk 


where is a proper numerical representation of the diffu¬ 
sive flux EE- The first partial update (step 4) in Sec. |3.3| ) 
for F then results to: 


F 2 = F + At* 


-(S t F) n hyp + ( StF ) 


"+2 

src 


cs -F" + CKtraAt*F^) — CKtraAf'i^ 
F CKtraAt* Jd 


1 + CKt ra A t* 1 + CKt ra At* 


(63) 


where At* = At/2. Hence, for CKtraAt* 3> 1 the flux density 
F consistently relaxes to Fr> within one (partial) time step 
without any numerical overshooting. 


3.7 Hydrodynamics 

The equations of hydrodynamics are integrated using the 
finite-volume high-resolution shock-capturing scheme devel¬ 


oped in Obergaulinger (2008). The scheme evolves the con¬ 


served hydrodynamic variables (p,pv,e t ), to which end a 
variety of procedures for spatial reconstruction (cf. Sec. 3.21 
and Riemann solvers (Lax-Friedrich, HLL, HLLC, HLLD, 
see, e.g., |Toro[l 997|) can be selected. Moreover, a magnetic- 
field solver and various models of viscosity and magnetic 
diffusivity are implemented. In the case of the model setup 
requiring the co-evolution of a set of different fluid species a 
simplified version of the “Consistent Multifluid Advection” 
scheme ( Plewa fe Miiller| |l999) is utilized. For more details 
about the non-radiative part of the code, we refer the reader 


to Obergaulinger (2008). 


4 TEST PROBLEMS 

In this section, we present a variety of problems to test the 
methods described in the previous sections. Several ideal¬ 
ized, non-microphysical tests in ID and 2D are performed 
which, although they are not directly related to typical sce¬ 
narios where neutrino transport plays a role, serve to assess 
the quality of the two-moment closure approximation and 
the coupling of the radiation moments to the velocity field 
and to the hydrodynamics part of the code. Subsequently, we 
present two one-dimensional test problems related to CCSNe 
in which we test the AEF scheme both for different closure 
prescriptions and against corresponding results from FLD 
and Boltzmann schemes. 

In order to avoid excessive repetitions we list some re¬ 
curring properties and parameters that various following 
tests are equipped with: 

• Regarding the numerical treatment, we employ a 5th- 
order MP reconstruction method and an HLL Riemann 
solver for both the radiation transport and the hydrody¬ 
namics part of the system. For the tests in ID we take a 
global CFL factor of CFL= 0.7 while for the 2D tests we 
set CFL= 0.5. 

• We apply boundary conditions (BCs) by fixing the values 
in the boundary (ghost) zones surrounding the compu¬ 
tational domain according to a given prescription. For a 
reflective boundary, e.g. at x = xo, we copy scalar quan¬ 
tities, e.g. E\ xq+Sx = E\ xq _ Sx , and apply the reflection 
operator to vectorial quantities, e.g. [F x , F v , F z )\ XQ+gx = 
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Velocity Fields 
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Spectral Differences 



Figure 1. Doppler shift of free-streaming radiation. Panel (a) shows the three velocity fields as functions of the location x and panel (b) 
depicts the differences between the energy distributions in the two frames with ft = 0,/3 max . The lines denote the fully relativistic, 
analytic solutions, while the symbols show the AEF results. 


(— F x , F v , F z )\ xq 5x . For a non-reflective outflow bound¬ 
ary, e.g. at * = xo we employ the usual Oth-order extrap¬ 
olation for all quantities, e.g. we set E(x o + Sx) = E(x o). 
In all of the subsequent tests in which dimensionless equa¬ 
tions and quantities are employed the speed of light is set 
to c = 1 and for the velocity the symbol P is used. 

Except for the tests in Secs. |4. l~3l |4.3. l| and |4.3(2| we will 
exclusively use the Minerbo closure, which is expressed by 
the Eddington factor as in Eq. (28a I and the 3rd-moment 
factor as in Eq. 


given by 


E f 3 


1 (« 0 3 
s e se — 1 ’ 


where s = 



(64) 


which is valid for 0^/3^ 1. 

The differences between the spectra in the frames P = 
0 , (9max for both our numerical and the analytic solution are 
shown in panel (b) of Fig. [I] The Doppler shift is captured 
well by our scheme: The agreement with the analytic solu¬ 
tion converges for decreasing values of /3 max . For high veloc¬ 
ities, /3 max = 0.3, the 0[v/c) approximation leads to errors 
of about 10% with respect to the relativistic solution. 


4.1 One-dimensional idealized test problems 

4-1.1 Doppler shift of free-streaming radiation 

The energy-coupling scheme describing the Doppler shift 
(cf. Sec. |3.5.3| ) can immediately be tested by comparing the 
spectra of a free-streaming radiation field in two different 
frames with non-vanishing relative velocity. We adopt the 
basic setup from |Vaytet et ah] (2011), but we use dimension¬ 
less units and take higher values for the maximum velocities 
/3 max £ {0.01,0.1,0.3}. The Cartesian spatial domain covers 
x € [0,10] and is resolved by 100 equidistant grid points, 
while the energy space is discretized between e G [0, 50] us¬ 
ing a logarithmic grid with 40 bins and a bin-enlargement 
factor of Ae^+i/Ae^ = 1.1. At x = 0 we inject a beam of 
radiation by fixing the radiation quantities in the boundary 
zones according to E(x = 0) = e 3 /(e e — 1) = F{x = 0). The 
boundary at x = 10 is set to outflow. The radiation field 
traverses velocity fields with the shape of smoothed step- 
functions as shown in panel (a) of Fig. [l] Within regions 
where P > 0 the redshifted (i.e. Lorentz-boosted) spectrum 
of the comoving-frame specific energy density is analytically 


4-1.2 Differentially expanding isothermal atmosphere 


To test the algebraically closed two-moment transport in 
combination with frame-dependent effects and (idealized) 
radiation-matter interactions we examine a scenario that 
was also investigated by Mihalas ( 1980| ) (in full relativity) 
and by Rampp & Janka (20021 (in 0(v/c)), both using ac¬ 
curate Boltzmann techniques. The scenario includes an ex¬ 
panding, isothermal atmosphere that expands with velocity 
beta as function of radius r as 


P{r) = Pu 


(65) 


and which exhibits an absorption opacity K a that varies in 
r and (photon) energy t as 


K a (r, e) = 


!§p e -(e-eo) /A + 


^ ( 


1 — e 


-( £ -e 0 ) 2 /A 2 


e ^ eo 


r , e > eo ■ 

( 66 ) 

That is, for fixed radius r, the opacity is a smoothed 
step-function in energy space with the transition at en¬ 
ergy eo from a low opacity to a 10 times higher opac¬ 
ity with transition width A. The model parameters are 
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Frequency Integrated Energy Density Spectral Distribution 




Figure 2. Differentially expanding isothermal atmosphere. In panel (a) the frequency-integrated energy densities, normalized by Eo = 
E(r = 0), as function of radius are shown for different maximum velocities (3 max . The symbols denote the reference solution computed 
by |Mihalas| ( |1980| , where the points are read off their Fig. 5. In panel (b) we show for the same models the spectral distributions of the 
normalized energy densities at the two radii r = 5.5 and 11, at which the optical depth for photons with energy e < eo is ~ 1 and ~ 0, 
respectively. The thin solid line displays the equilibrium distribution F eq (e). 


{ f 'min ) T max 5 eo,A,a} = {1,11,3,0.2,10.9989} and the max¬ 
imum velocity /3 max is varied between /3 max € {0,0.1, 0.3}. 
We use dimensionless units in which the temperature T = 1 
such that the photon equilibrium energy density is given by 
_E eq (e) = e 3 /(e e — 1). We set up a uniform radial grid of 400 
zones to cover a region of r £ [0.1,15]. The additional (with 
respect to the reference calculations) region [11,15], wherein 
the opacities are set to zero, merely serves as a transition 
zone for the radiation field to reach near free-streaming con¬ 
ditions in order to avoid unphysical feedback from the outer 
radial boundary, at which outflow BCs are employed. At 
r = 0.1, a reflective BC is applied. The energy grid is com¬ 
posed of 50 equidistant bins within e £ [0,12], 


well by the AEF method. We notice an increasing difference 
between the 0{v/c ) results and the fully relativistic results 
for increasing /3 ma x. However, for /3 ma x = 0.3 the maximum 
error in the integral energy density, cf. panel (a) of Fig. [2j 
is still less than ~ 5 %. 


4-1-3 Supercritical radiative shock 

Successively increasing the degrees of freedom taken into 
account, we now turn to a classical RHD problem to test 
the accurate coupling between transport and hydrodynam¬ 
ics, the radiative shock tube. Having been the subject of 


numerous investigations, both analytically (e.g Zeldovich & 


sity as function of radius is shown in panel (a), while the 
radiation spectrum at two representative radii is shown in 
panel (b). A remarkable fact is that the case with no expan¬ 
sion is already reproduced well with an accuracy of < 1% 
by the approximate AEF scheme. By switching to P > 0 
we introduce the following effects: Due to the expansion the 
comoving-frame energy- and flux densities of photons cre¬ 
ated deeper within the atmosphere decreases on their way 
to the surface, as can be seen by monotonic decline of the en¬ 
ergy densities with increasing /3 ma x up to r ~ 10. However, 
this trend is competed by the effect that the overall frac¬ 
tion of photons originally created in the high-opacity band 
is lifted with increasing velocity (and radius), since a frac¬ 
tion of photons are redshifted from e > eo to e < eo on their 
way to the surface. Hence, the opacity jump is effectively 
redshifted by the expansion, leading to higher integral val¬ 
ues of the energy density at r > 10. Both effects are captured 


For comparison, we show similar plots as in Rampp &| 

Raizer 1966 Mihalas &; Mihalas|1984) and numerically (e.g. 

Janka (2002), see Fig. |2| The total radiation energy den- 

Ensman 1994 Sincell et al. 19991, radiative shock tubes re- 


peatedly served as test problems for the development of new 


RHD codes as, e.g., in Turner & Stone 

(2001); Hayes & Nor- 

man (2003); Gonzalez et al. (20071; 

Vaytet et al. (2011); 

Jiang et al. (20121; Skinner & Ostriker (|2013). 


Since the detailed physical description of radiative 
shocks is out of the scope of our presentation, we only briefly 
summarize their essential properties here. In contrast to 
purely hydrodynamic shocks radiative shocks allow for en¬ 
ergy transfer between the gas and radiation, effectively in¬ 
troducing cooling of the post-shock and heating of the pre¬ 
shock material. Depending on the shock velocity, heating of 
upstream material in front of the shock - this region is called 
radiative precursor - can become so efficient that the pre¬ 
shock temperature adapts to the post-shock temperature, 
in which case the shock is called a supercritical radiative 
shock. In this case, both the up- and downstream material 
is in radiative equilibrium close to the shock and separated 
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Radiative Shock 



0 2 4 6 8 10 

x [ 10 10 cm] 


Figure 3. Supercritical radiative shock. We show the gas and 
radiation temperatures at three different times t £ {4 X 10 3 , 7.5 X 
10 3 ,1.3 X 10 4 } s. Each curve is plotted in the frame in which the 
shock crosses x = 0 at t = 0 with a velocity »,=2x 10® cm „s~ 1 . 
The square symbols denote points from the reference calculation 
by | Vaytet et al.| (2011), which have been read off their Fig. 9. 


by a sharp non-equilibrium temperature spike (“Zeldovich 
spike”), which is roughly as wide as the local mean-free path 
of radiation. 

We initialize our model of a supercritical radiative shock 
using the same setup as Vaytet et al. (2011), who also 
employed a spectral AEF scheme fairly similar to ours. 
Moreover, to facilitate the comparability with [Vaytet et ah] 
(2011) we apply the Mi instead of the Minerbo closure for 
this test. The only difference to the setup of [Vaytet et al.j 
(2011) is that we ignore here the 3rd-moment terms in 
the evolution equation for the 1st moment (i.e. the (IV)- 
terms in Eq. ©)■ However, these terms vanish in the 
energy-integrated form of the RHD equations and the fluid 
quantities, such as the temperature, should not be sig¬ 
nificantly influenced by this measure, particularly in view 
of the relatively low velocities at hand. A Cartesian box 
of length 10 11 cm is discretized by a uniform grid of 500 
cells and initially homogeneously filled with gas of density 
p = 7.78 x 10 -10 gcm -3 , temperature T = 101< and grey 
absorption opacity K a = 3.1 x 10~ 10 cm _1 . Initially, radi¬ 
ation is everywhere in equilibrium with the gas, such that 
E = E eq = OradT 4 , where a ra d « 7.57 x 10~ 15 ergcm~ 3 K -4 
is the radiation constant. The gas pressure is computed as 
P s = (7gas — l)e; = pk B T/mB with 7 gas = 1.4. We take 
the frame of the shock moving with v s = 2 x 10 6 cms _1 
relative to its preceding medium as our simulated inertial 
frame. To this end, we let all matter in the computational 
domain initially move with velocity v = — v s . The shock is 
induced by using a reflective boundary at x = 0 and it is 
maintained by feeding new material with the original prop¬ 
erties at x = 10 11 cm into the computational domain. We 


discretize the frequency space with N e = 8 evenly spaced 
bins between e £ [0, 8 x 10 14 ] Hz. 

The results are shown in Fig. [3] in form of the distri¬ 
butions of the gas temperature T gas and the radiation tem¬ 
perature, the latter being defined as T ra d = (E/ar^d) 1 ^ 4 ■ Us¬ 
ing an essentially similar physical evolution model as|Vayte~t| 


et al. (2011), the results we obtain with our quasi-explicit 


numerical method are in good agreement with the outcome 
of their implicit radiation solver which shows that the cou¬ 
pling between the radiative and hydrodynamics systems is 
numerically robust and produces accurate results. 


4.2 Two-dimensional idealized test problems 


4-2.1 Shadow casting problem 

We begin our presentation of 2D test problems with a rather 
qualitative test that puts one of the generic advantages of 
a two-moment scheme into focus, namely the ability of an 
opaque object to generate a shadow when being illuminated 
by radiation. In a two-moment scheme, the flux-density vec¬ 
tor is an evolved quantity, and the pressure tensor is gen¬ 
erally non-isotropic. In contrast, in an FLD scheme (see 
Sec. |2.4.1[ ), the flux direction is determined by the gradi¬ 
ent of the scalar energy density, which corresponds to an 
isotropic radiation pressure. This leads to the unphysical ef¬ 
fect that in the free-streaming regime, / —> 1, anisotropic 
features of the radiation field, such as shadows, cannot be 
maintained and are quickly washed out. 

As it has likewise been done before for various other ra¬ 
diative transfer/transport codes (see e.g. Audit et al.||2002 


Hayes fe Norman|2003| |Iliev et al.|2006| |Skinner fe Ostriker| 

2013| and references therein) we set up a purely absorb¬ 
ing gas cloud which is exposed to near free-streaming radia¬ 
tion to test the ability of radiation to cast a shadow behind 
the gas cloud. Specifically, in a Cartesian domain with x = 
(x, y) £ [0,15] x [—5, 5] and resolved by N x x N y = 300 x 200 
cells, we define one region in which the radiation field is gen¬ 
erated, the circular region S centered around xs = (3, 0) 
with radius rs = 3/2, and we define another circular region 
A centered around xa = (11,0) with radius ta = 2 to be 
the purely absorbing cloud. The absorption opacity « a and 
equilibrium energy density _E eq are defined as follows: 

, x £ S 

, x £ A (67a) 

, else , 

(67b) 

The model is initialized with vanishing flux densities and a 
homogeneous distribution of negligibly small energy densi¬ 
ties. 

From the numerical point of view, the present objective 
is to test the correct implementation of the multidimensional 
hyperbolic part of the radiation moment equations, partic¬ 
ularly of the angular dependence of the signal speeds in the 



Riemann solver, cf. Eqs. (471. The signal speeds determine 


the numerical fluxes between grid cells, cf. Eq. (481, and 


close to free-streaming conditions both the signal speeds as 
well as the intercell fluxes should be strongly suppressed or¬ 
thogonal to the direction of the radiation flux. 
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Figure 5. Static and dynamic diffusion. We show profiles of the energy densities along the slice y = 1 after the simulation time t — to = 5 
of the Gaussian pulse centered at xq = (1,1). In panel (a) the results for the cases with spatially constant opacity are plotted. The models 
differ in spatial resolution, and therefore in the stiffness parameters V, and in the uniform velocity field parallel to the x-direction with 
magnitude 0. The results for the cases with 0 > 0 were transformed to the 0 = 0 frame and shifted to their initial position at t = 0. For 
a better readability, only every second zone value is plotted for the high-resolution models. The results shown in panel (b) are produced 
by the same initial Gaussian pulse but employing an opacity that declines with distance to the center of the pulse. The vertical lines 
denote the locations where the stiffness parameter of the low-resolution models is "P = 1. In the high-resolution models, the stiffness 
parameter is V < 1 everywhere. 


The three snapshots in Fig. [4] show for three consec¬ 
utive times contours of the isotropic luminosity L emitted 
by the source, given in this two-dimensional geometry by 
L = 2nr c \F\ with r c = \x — ats|. One can see that a clearly 
obscured region behind the gas cloud emerges. The lumi¬ 
nosity behind the gas cloud is not an ideal step-function in 
vertical direction but it changes rather continuously within 
a fan of opening angle « 20° — 30°. The reasons for this 
are, first, that radiation is not emitted from a point-like but 
a spatially extended source, causing the flux-factor to be 
|-F|/.E ~ 0.98 < 1 at r c = 8, and second, that the gas cloud 
is not perfectly absorbing but has a finite value of K a , which 
allows a small fraction of radiation to pass through the gas 
cloud near its edges. Altogether our code performs well in 
this test, the development and propagation of the multidi¬ 
mensional radiation field and its particular feature to cast a 
shadow are consistently captured. 


4-2.2 Static and dynamic diffusion 


A standard test for radiation codes allowing for the treat¬ 


ment of optically thick regions (e.g. Gonzalez et al. 2007 


Swesty & Myra 20091 is the scenario of an initially con¬ 


centrated bulge of radiation diffusing into its environment. 
Being conceptually based on the diffusion limit, an FLD 
scheme usually performs well in this test, as long as the 
medium is sufficiently opaque. On the other hand, in a two- 
moment AEF scheme the diffusion equation only results 
in the parabolic limit of the otherwise hyperbolic equa¬ 
tions. Therefore, it has to be checked that the numerical 


method chosen to solve the AEF scheme consistently de¬ 
scribes the parabolic diffusion limit and that the transition 
to the hyperbolic regime proceeds in a numerically stable 
and accurate fashion (cf. Secs. 3.4 and 3.6 1 . Concerning the 
last point, we explicitly want to ensure that no spurious, 
resolution-dependent features result from modifying the nu¬ 
merical fluxes whenever the local stiffness parameter exceeds 
unity, cf. Eq. (501. Finally, in this setup we also want to test 
the ability of the algorithm to accurately describe dynamic 
diffusion, i.e. diffusion out of a moving medium. 

We perform a set of calculations in a Cartesian box 
given by x = (x,y) € [0,3] x [0,2]. All configurations are ini¬ 
tialized at a fiducial time to = 5 with the following Gaussian 
pulse of radiation energy density and corresponding diffusive 
flux density centered around Xq = (1,1): 


E(x,to) = Eq exp 


- 4 ^° |2 } , F{x,t 0 ) = -D 0 VE, 


( 68 ) 

where Eq = 1 and Do = 3 x 10 3 . 

In the first configuration we define a spatially constant 
diffusion coefficient, D = (3k s ) _1 = Do, which allows us 
to compare the numerical results with an analytic solution, 
given by: 


E(x,t) = Eo ——— exp 
V ; t 0 +t 


\x— Xp \ 2 \ 

4Z?o(to + 1) J 


(69) 


For this configuration we switch between two resolutions 
{N x ,N y } = {450,300} and {150,100} corresponding to 
which the stiffness parameters, V = k b Ax, are lower and 
greater than 1, respectively. Additionally, for both resolu- 
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K? 


Velocity Field 


Rel. Difference of Flux at r = 10 
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(b) 


polar angle [7t] 


Figure 6. Radiation traversing variable velocity fields. In panel (a) we show the absolute value (color coded) of the velocity and its 
polar component (arrows) for the high-velocity case with /3 max = 0.1. In panel (b) we plot the relative difference of the energy-integrated 
radial flux density, F r , with respect to the case of vanishing velocity, F x . measured at radius r = 10 as function of polar angle. 


tions we vary between a vanishing and a non-vanishing spa¬ 
tially constant velocity in ^-direction, i.e. /3 = 0 and 0.1, 
respectively. For comparison with the /3 = 0 models, we plot 
the radiation energy densities of the /? > 0 models in the 
coordinate frame of the /3 = 0 models. 

In panel (a) of Fig. [5] we see that after a simulation time 
of t — to = 5, when the maximum value of E has reached 
about half of its initial value Eq = 1, all models agree well 
with the analytic solution and we cannot identify unphysical 
numerical features in any of the different cases. 

We consider a second configuration with the same initial 
conditions, cf. Eqs. ( 681 , to test the correct transition from 
a stiff (V > 1) to a non-stiff (V < 1) regime. To this end, 
instead of a constant opacity we now use a variable opacity 
that declines with distance from the center as 


k s (x) 


1 

3Do 


exp 


\x 



(70) 


on a length scale of S = 0.4. We display the results at 
the simulation time t — to = 5 in panel (b) of Fig. [5] We 
again compare four cases with two resolutions, {N x ,N y } = 
{450,300} and {225,150}, and two homogeneous velocities 
in ^-direction ft = 0 and 0.1. For the models with /3 > 0, 
the opacity profile is advected with the fluid, i.e. the opac¬ 
ity Ks(x,t) = k s (x— fit) is employed. In the high-resolution 
cases, which serve as references for the low-resolution cases, 
the stiffness parameter is V < 1 everywhere, while in the 
low-resolution cases the stiffness parameter crosses V = 1 
at the locations indicated by the dotted lines. We observe 
no numerical artifacts near the transition where V = 1 in the 
low-resolution cases, which indicates that the modification 
of the numerical fluxes, Eq. (501, works accurately. 


4-2.3 Radiation traversing variable velocity fields 

Since the quantities E, F evolved in our two-moment formu¬ 
lation are defined in the comoving frame, they are subject to 
variations whenever radiation crosses regions of variable ve¬ 
locity even without any interactions present. The net impact 
on the radiation properties after passing such regions and re¬ 
turning back into the original frame would vanish in an exact 
calculation. In practice, however, we encounter two limita¬ 
tions that spoil this feature to be fulhlled precisely: First, 
our underlying scheme for the radiation moments neglects all 
contributions of order 0(v 2 /c 2 ) in both evolution equations, 
which results in a loss of the property that a transformation 
from one frame to another is exactly reversible - instead 
such a transformation generates errors of the disregarded 
order 0(v 2 /c 2 ). The second reason is that we do not solve 
the evolution equations exactly but only numerically, i.e. all 
our solutions are beset with truncation errors depending on 
the spatial and temporal resolution and on the numerical 
algorithm. 

In order to obtain a qualitative impression of how 
strongly both aforementioned effects can impact the solu¬ 
tion, we set up an arbitarily shaped velocity field and let 
it be traversed by a spherically expanding radiation field. 
The radiation field is induced at the boundary at radius 
r = 2 with an energy density of E(e) = e 3 /(e £ — 1), a ra¬ 
dial flux density of F r (e) = E(e)/2, and vanishing non-radial 
flux components, Fg = = 0. The energy grid consists of 

N c = 10 bins logarithmically distributed between e € [0, 30] 
with an enlargement factor Ae^+i/Ae^ = 1.3. The velocity 
field in the polar plane, v po i, represents an eddy with radius 
di = 1 circulating around its center at xo = (5, 5)/v2, while 
the toroidal velocity field v t0 r has the same absolute mag¬ 
nitude as the poloidal field but points into the (^-direction. 
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Figure 4. Shadow casting problem. The 2D luminosity at the 
three indicated times t is plotted. The dashed line indicates the 
boundary of the absorbing gas cloud. 


Explicitly, we set 


_ P(x) _ f}(x) 

Up oi — epo1 1 VtoT ~ 2 


(71) 


with 


/3(x) = /3 max exp{-(|a;- xq\ - rfi) 2 /rf|} , (72) 


where d 2 = 0.4, is the unit vector in ^-direction, and 
e po i is the unit vector perpendicular to both x — xo and 
and signed such to point in clockwise direction in the R — z 
plane (see panel (a) of Fig.|6|. We vary the maximum value 
of the velocity between /3 max £ {10 —3 , 10 —2 , 10” 1 } and we 
use the two spatial resolutions N r = Ng £ {50,100} between 


r £ [2,10] and 6 £ [0, tt/ 2] . At r = 10 we employ an outflow 
BC. 

In panel (b) of Fig. [6] we compare the energy-integrated 
fluxes F r obtained for each velocity field and resolution 
with the corresponding value F° resulting for a vanishing 
velocity field. If the two shortcomings mentioned in the 
beginning of this section were absent, both fluxes would 
be exactly equal, i.e. (F r — F°)/Fr = 0. Instead, in our 
numerical calculation we receive relative errors of up to 
{F r -F?)/F° r | max ~{4x 10 -2 ,4 x 10- 4 ,1 x 1(T 5 } for 
/3max = {10 -3 ,10 -2 ,10~ 1 }. The outcome that the lines lie 
much closer together for different resolutions than for differ¬ 
ent 13 max shows that the low resolution already sufficiently 
resolves the solution corresponding to the underlying mo¬ 
ment equations, at least for /? max > 10~ 3 . It can further 
be observed that the leading-order error term representing 
missing components compared to the fully relativistic formu¬ 
lation is not 0{v/c), as can be inferred from the tendency of 
relative differences to roughly decrease by two orders of mag¬ 
nitude for a one-order reduction of /3m ax- Hence, we deduce 
that in our implementation no significant contributions of 
order 0(v/c ) are missing or are erroneously present since in 
any other case we would have found an error scaling linearly 

with /fmax- 


4.3 One-dimensional problems including 
microphysics 

While the previous test problems are based on rather ide¬ 
alized setups, the two remaining one-dimensional test prob¬ 
lems specifically focus on neutrino transport in CCSNe. In 
particular, the tests should address the question how the 
AEF scheme performs compared to the present standard 
methods, the FLD and the Boltzmann-type solvers. To this 
end, in the first test in Sec. |4.3.T] we keep the hydrodynamic 
background - consisting of a typical proto-NS configuration 
- fixed and we only compare the stationary radiation field 
resulting from the three different aforementioned types of 
methods. In the second test in Sec. |4.3.2| we then compare 
the results of a fully dynamic CCSN simulation with two 
similar calculations performed with well-known Boltzmann- 
type neutrino-hydrodynamics codes. 


4-3.1 Neutrino radiation field of a static proto-neutron 
star 

To compare the different neutrino transport schemes AEF, 
FLD and Boltzmann with each other in the CCSN context in 
a manner that is independent of the hydrodynamics part of 
the numerical method, it is instructive to evolve only the ra¬ 
diation field in a proto-NS background that is held constant 
during the evolution. As background configurations, we take 
two profiles of hydrodynamic data obtained from two differ¬ 
ent simulations performed with the VERTEX code. The hy¬ 
drodynamic data for our first model (called “pre-explosion 
model” hereafter) is represented by a snapshot taken at time 
150 ms after bounce in the accretion phase of model “N13”, 
which was investigated in the course of |Liebendorfer et al.| 
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Pre-Explosion Model Post-Explosion Model 



Figure 7. Neutrino radiation field of a static proto-neutron star. We compare for two models (see top title) between different transport 
schemes (see line labels in the second panel from the top). In the panels are displayed from top to bottom the properties (density p , 
temperature T and electron fraction Y e ) of the hydrodynamic background, the mean flux-factor / = F/(cE), luminosity L = Airr^F 
and rms-energy (e) rm s = yjf eE(e)de/ f N(e)de of electron neutrinos, and the source terms Qe/^B? Qn /P for the gas energy density 
(cf. Eq. ( |l3a[ )) and electron-number density (cf. Eq. < |l3c[ |), respectively. Note that in cases when the dotted green line is invisible it is 
superimposed by the solid green line. 
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( 2005|p^ The snapshot for our second model (called “post¬ 


explosion model” herea fter) is taken from the model “Sr” in 
Hiidepohl et ah! (2010|j^ in the neutrino-driven wind phase 


at time 300 ms after core bounce. In both models we adopt 
the spatial and energy grids and the profiles of the density 
p, temperature T and electron fraction Y e from the reference 
calculations, while the remaining quantities needed for the 
opacities have been obtained by applying the EOS of |Lat-| 


timer & Swesty (1991). Note that since the data in the pre¬ 


explosion model is partially rather noisy we smoothed the 
original profiles of p°, i\ V)° as given by the reference calcu¬ 
lation by making the replacement u, = (uP_! +2 uf -|- ■up ) _ 1 )/4 
for u € {p, T, Ye} at each radial grid point i in the pre¬ 
explosion model. In both models, we ignore frame-dependent 
effects by setting the velocity v = 0 everywhere. Regarding 
the neutrino interactions, we only take into account the nu¬ 
cleonic /3-processes (n-|-e + p+h e andp+e - *=; n+u e ) and 
isoenergetic scattering of (anti-)neutrinos off free nucleons as 
described in Bruennj (1985). Only eletron-type neutrinos are 
evolved. 

For comparison of the AEF scheme with the two re¬ 
maining schemes, we additionally implemented an FLD 
scheme and a Boltzmann solver. The FLD scheme was 
obtained by modifying our existing time-dependent AEF 
scheme while for the Boltzmann-type solution we imple¬ 
mented a separate algorithm in which the time-independent 
two-moment system is solved with the closure being pro¬ 
vided by a tangent-ray scheme. The details about both 
methods can be found in Appendix [C| 

For each of the two background models eight calcula¬ 
tions using different methods or closures were performed: 
The reference solution is represented by the calculation 
conducted with the Boltzmann-type tangent-ray variable- 
Eddington-factor scheme (TR-VEF), since this method 
should yield the most accurate solution. In four additional 
simulations AEF schemes together with the Eddington fac¬ 
tors in Eqs. (28) were employed, while in the remaining 
three simulations FLD schemes together with the flux- 
limiters in Eqs. (251 were applied. For computing the flux- 
limiter ABruenn, cf. Eq. (25bI, the (energy-dependent) neu- 
trinosphere radii r„, defined by 


r oo 

r v [r v , e)= Ktra(r', e)dr' = 1, 
J r ,, 


(73) 


were used. 

In Fig. [TJ we compare for the different calculation meth¬ 
ods the radial profiles of the mean flux-factor, luminosity, 
and rms-energy of electron neutrinos, as well as the hydro- 
dynamic source terms corresponding to heating/cooling and 
(de-)leptonization resulting from the interaction with both 
electron neutrinos. We observe the following properties: (1) 
Concerning the luminosities and flux-factors, the accuracy 
of the AEF schemes is throughout slightly higher than that 
of the FLD schemes. While the FLD luminosities and flux 
factors exhibit errors up to ~ 10%, the corresponding AEF 
quantities are throughout accurate to less than ~ 3 %. (2) 


The rms-energies are accurate to <( 1 % in all schemes. (3) 
The accuracy to which the hydrodynamic source terms are 
reproduced is comparable for AEF and FLD in the pre¬ 
explosion model, with the heating rates in the gain region 
being somewhat underestimated for FLD and overestimated 
for AEF. In contrast, in the post-explosion model the AEF 
schemes perform clearly better than the FLD schemes. The 
main difference between both background models is that 
the semi-transparent region is more extended in the pre¬ 
explosion model. (4) The overall results appear to be more 
sensitive to different flux-limiters when using FLD than to 
different Eddington factors when using AEF. 

The results of this comparison indicate that AEF meth¬ 
ods can perform at least equal to, if not better than, FLD 
schemes in a ID proto-NS environment. However, since we 
have only investigated two stationary snapshots we cannot 
be sure about the universality of the observed levels of ac¬ 
curacy with respect to time- and model variations. For this 
reason we avoid at this point to speculate about the hydro- 
dynamic impact of the observed differences in the heating- 
and deleptonization-rates for the various methods. Never¬ 
theless, the comparison test discussed in the next section 
demonstrates that AEF schemes can in fact compete with 
Boltzmann solvers even in the fully dynamic case. 


4-3.2 Fully dynamic collapse and post-bounce evolution of 
a 13 Mg progenitor 

Our final and most comprehensive test comprises the self- 
consistent collapse and subsequent post-bounce evolution 
of the core of a star with approximate main-sequence 


mass of 13 Mg (Nomoto & Mashimoto 1988). This model 


was already investigated in the comparison in |Liebend orfcr 
et al. ( 2005|) (labeled “N13” there) between the neutrino- 


hydrodynamics codes AGILE-BOLTZTRAN (Liebendorfer 


et al. 

20041 

Janka 

2002) 


2004) and VERTEX-PROMETHEUS (Rampp & 


2002). The comparison revealed a good agreement 


between the two Boltzmann-type codes, which allows for a 
straightforward assessment of the accuracy of our method 
by dir ect comparison with the r esults of Liebend orfer et aT7| 
12005 Note in passing that O’Connor (20141 recently 
conducted a similar comparison but using the AEF-type 
code GRID and inspecting model “G15” fro m |Liebendorfer| 
et al. (20051. The main finding of O’Connor ( 2014| ), namely 
the very good agreement of the AEF scheme with the Boltz¬ 
mann codes, is in consensus with ours. Finally, using ad¬ 
ditional, slightly modified simulation setups we will check 
the robustness of our results with regard to the numerical 
scheme and the choice of the closure prescription. 


4.3.2.1 Model setup We start the spherically symmet¬ 
ric evolution at the onset of collapse when the core has 
a central density of p c = 3.16 x 10 10 g cm -3 and we fol¬ 
low the system up to a post-bounce time of t = 300 ms. 
This test involves all types of neutrino interactions listed 
in Sec. |2.3| Thus, except for our omission of the processes 


13 See Sec. |4.3.2| for more details about this model, which also 
served as reference model in the fully dynamic CCSN simulation 
presented in that section. 

14 The data was provided via 

http://www.mpa-garching.mpg.de/ccsnarchive/ . 


15 The results discussed in iLiebendorfer et all 
(2005| and used here as reference data are pub- 
lically available under http://iopscience.iop.org/0004- 
637X/620/2/840/fulltext/datafiles.tar.gz 
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Figure 8. Collapse dynamics in the core-collapse test. All black, orange and green lines display results obtained with our code, VERTEX 
and AGILE-BOLTZTRAN, respectively. The panels (a) and (b) show the evolution of the central electron (solid lines) and lepton fraction 
(dashed lines) and the central entropy, respectively, as functions of the central density during collapse. The increase of the central entropy 
for central densities p c > 10 12 gem -3 in the AGILE-BOLTZTRAN run is a numerical artifact which has no impact on the subsequent 
physical evolution (see |Liebendorfer et al. 2005| for more details). The panels (c) and (d) display profiles of the electron fraction and the 
entropy as functions of the enclosed mass at times when the core reaches the central densities p c shown in the legend in panel (c). 


i/ e P e *=> e~e + P^| we empl oy the same neutrino microphysics 
as in |Liebendorfer et al. ( |2005| ). Above a threshold density 
of Pls = 6 x 10' g cm" J the stellar gas is modeled by the 
nuclear EOS of Lattimer & Swesty ( |1991 1, with a compress¬ 
ibility modulus of K = 180 MeV. Below this threshold, we 
use a low-density EOS that includes photons, arbitrarily rel¬ 
ativistic and arbitrarily degenerate electrons and positrons, 
and a non-relativistic Boltzmann gas of baryons. For p > pLS 
the nuclear composition is in nuclear statistical equilibrium, 
which only depends on p,T and Y e . For p < pls we assume 
the baryonic composition to be given by pure 28 Si for tem- 


1,1 However, these reactions are subdominant in both the collapse 
and post-bounce evolution by at least an order of magnitude com¬ 
pared to the dominant reactions, see e.g. IBuras et al. 12006 


peratures T < 0.44 MeV, or by pure 66 Ni for temperatures 
T > 0.44 MeV. 


We employ a Eulerian radial grid with N r = 384 zones 
distributed logarithmically between the origin and an outer 
radius of ~ 7, 000 km with a minimum grid width of 200 m. 
The neutrino energy space is discretized into 21 energy bins 
that are roughly logarithmically distributed between 0 MeV 
and 400 MeV. For the initialization of the model we take 
the profiles of velocity, density, electron fraction and entropy 
from the corresponding AGILE-BOLTZTRAN run. For the 
time integration of the source terms we ignore the implicit 
time dependence of the hydrodynamic quantities. Specih- 
call y, we use the integration modes c) and b) described in 
Sec. 


3.6 


for densities lower and higher than 10 11 gem 3 , 


spectively. With the CFL factor set to 0.7, our whole simu¬ 
lation required the calculation of about 700, 000 time steps. 
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Hydro-Profiles, t - 10 ms Neutrino-Profiles, t - 10 ms 



Hydro-Profiles, t = 150 ms Neutrino-Profiles, t - 150 ms 



Figure 9. Radial profiles in the post-bounce phase of the core-collapse test. All black, orange and green lines display results obtained 
with our code, VERTEX and AGILE-BOLTZTRAN, respectively. In panel (a) we show the density and velocity (top), and the entropy 
per baryon and electron fraction (bottom), while in panel (b) we plot the luminosities (top) and rms-energies (bottom) at a post-bounce 
time of t = 10 ms. In panels (c),(d) the same respective quantities as in panels (a),(b) are plotted, but at a post-bounce time of t = 150ms. 


4.3.2.2 Collapse dynamics The core collapses in a 
time of t co ii ~ 95 ms to reach a maximum density of 
Pm ax = 3.6 x 10 14 g cm~ 3 . We show the evolution of the 
central values of the electron fraction Y e , lepton fraction 
Viep = (n e - — n e + + N Ut , — Nv e ) /ub and entropy per baryon 
s as functions of the central density p c during collapse in 
the top panels and the profiles of Y e and s as functions 
of the enclosed mass coordinate M enc (r) = 4-7T fj pf 2 dr in 
the bottom panels of Fig. [8] respectively. Starting at central 
densities of about p c ~ lCr 1 g cm -3 , inelastic electron scat¬ 
tering reduces the mean energy of the neutrinos escaping 
from the core, leading to an accelerated deleptonization and 


to an increase of the central entropy of the gas. Neutrino 
trapping sets in at a central density of p c sw 2 x 10 12 g cm -3 , 
above which the central values of Y) ep and s remain roughly 
constant until core bounce. In all variables, we find a nearly 
perfect agreement of our results with both reference solu¬ 
tions. 

4.3.2.3 Post-bounce evolution Once the core reaches 
nuclear densities, a shock wave is formed at a mass coor¬ 
dinate of M s h « 0.67 Mg, consistent with the reference re¬ 
sults. We show the radial profiles of various quantities at 
two different post-bounce times in Fig.[9]and the time depen- 
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Shock Radii Luminosities / Mean Energies 




Figure 10. Post-bounce evolution in the core-collapse test. All black, orange and green lines display results obtained with our code, 
VERTEX and AGILE-BOLTZTRAN, respectively. In panel (a) we show the shock radii as functions of the post-bounce time. The top and 
bottom of panel (b) depict the luminosities and mean energies, respectively, of electron neutrinos (solid lines) and electron antineutrinos 
(dashed lines). Note that for the luminosities the scaling on the left y-axis applies for t < 50ms, while the scaling on the right y-axis 
applies for t > 50 ms. 


dence of the shock radius as well as the neutrino luminosities 
and rms-energies in Fig. m All quantities are defined as in 
Liebendorfer et al. ( |2005 1. About m 5 ms after bounce, the 
prompt shock reaches a radius of r s h oc k ~ 100 km, where it 
stalls for almost 10 ms. At this point the slightly stronger 
prompt shock in VERTEX leads to somewhat larger shock 
radii in the VERTEX runs compared to our simulation and 
to the AGILE-BOLTZTRAN run. Subsequently, neutrino 
heating increases the entropy of the matter in the gain re¬ 
gion behind the shock and causes the shock to propagate 
further outward. However, the energy deposition by neutri¬ 
nos is not sufficient to drive an explosion and we observe a 
continuous decrease of the shock radius after it has reached 
its maximum value of r s hock ~ 250 km at t ~ 120 ms. The 


shock trajectory r s hock(t) (cf. panel (a) of Fig. 101 of our 
simulation is in good but not perfect agreement with the 
reference results which, however, are also not exactly con¬ 
sistent with each other. The quantitative differences could 
to some degree stem from slight differences of the thermo¬ 
dynamic treatment of the low-density regime p < pLS to 


which Liebendorfer et al. (2005 I already attributed the dis¬ 


crepancies of r s hock(t) between the two reference models: 
Differences in the treatment of nuclear burning in the low- 
density regime lead to different entropies of the material that 
falls into the shock and thereby to different post-shock en¬ 
tropies. Besides further potentially significant differences in 
the hydrodynamic treatment of the three codes, the small 
underestimation of r a hock(l) in our AEF run at late times 
could also be the result of the approximate closure. However, 
even though the situation remains unclear, we do not con¬ 
sider this behavior to be very significant, given the fact that 
even both Boltzmann codes exhibit differences in r s hock(t) 
of roughly the same size. 


4.3.2.4 Neutrino emission At bounce, the core emits 
a short, intense burst of electron neutrinos. The neutrino 
burst in our simulation exhibits a higher maximum lumi¬ 
nosity of L max « 5.3 x 10 53 ergs _1 than in the two ref¬ 
erence solutions (I/max « 4.14 x 10 53 ergs -1 and L max ~ 
4.55 x 10 53 ergs -1 for AGILE-BOLTZTRAN and VERTEX, 
respectively). However, the shorter duration of the burst in 
our model compensates for this, and the integrated energy 
emitted during the burst, Eburat = / 0 °'° 2s L Ve (t) dt, is al¬ 
most equal for all three codes: Eburst = 2.80, 2.80, and 
2.85 x 10 51 erg for our model, AGILE-BOLTZTRAN, and 
VERTEX, respectively. The reason for our neutrino burst 
being sharper may be found in the numerical treatment: 
Our explicit time integration (at least for the hyperbolic 
terms which describe the propagation of radiation) com¬ 
bined with a high-order spatial discretization is certainly 
less dissipative than both fully-implicit reference schemes, 
which facilitates the accurate evolution of narrow radiation 
peaks. Note that this explanation is supported by the fact 


that O’Connor (2014) obtain a similarly enhanced neutrino 


burst using their explicit AEF scheme; see their Fig. 6. 

After the end of the neutrino burst, the luminosity of 
electron neutrinos drops quickly while the luminosity of elec¬ 
tron antineutrinos increases, and subsequently neutrinos of 
both flavors are emitted at about equal luminosities. After 
the mean energies of both flavors peak at the neutrino burst 
they remain approximately constant with a slow trend to¬ 
wards higher values. Due to their larger interaction rates 
with matter, electron neutrinos decouple at lower densities 
and temperatures than electron antineutrinos, leading to a 
mean energy that is ~ 3 MeV below that of the antineutri¬ 
nos. 

In summary, most of our results lie well within the tol- 
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Shock Radii Luminosities / Mean Energies 




Figure 11. Same as Fig. |10| but for the runs “standard”, “altered numerics”, “Mi closure”, and “maximum entropy closure” described 
in text. The maxima of the neutrino bursts in these models are L max r; {5.3,5.3, 5.1, 5.4} X 10 s3 ergs -1 , respectively. 


erance region spanned by the results of both reference calcu¬ 
lations, such that no obvious deficiencies of the AEF method 
can be identified. 


4.3.2.5 Variations of the calculation method In this 
last study, we address various issues to check the robustness 
of the results presented above and therefore the reliability 
of the whole algorithm described in this paper. To this end, 
we perform additional test runs that are described below. In 
what follows, we denote the simulation presented above as 
the “standard” run. In Fig. [IT] we show similar quantities as 
in Fig.[l0]but we compare the “standard” run with the test 
runs described below. 

We set up a simulation identical to the “standard” 
run described above but with altered numerical specifica¬ 
tions, called “altered numerics”. The number of spatial grid 
points is increased to 576 (with an unchanged innermost 
grid width of 200 m) and the time step is reduced according 
to CFL=0.3. Moreover, the integration mode for the source 
terms a) is used instead of b) for densities p > 10 11 g cm -3 
(cf. Sec. 3.6 1 . As a final difference to the “standard” run, we 
additionally take into account all the velocity-/acceleration- 
dependent terms that have been dropped when deriving 
the moment equations, Eqs. 0 . from the transfer equa¬ 
tion, Eq. Q, in Sec. [5] For the numerical implementa¬ 
tion of these terms we adopt methods in close analogy to 
the existing ones, cf. Sec. [3] This model yields very simi¬ 
lar results compared to the standard case, as is exemplarily 
shown for selected functions of time in Fig. ED This sim¬ 
ilarity between the two models proves the robustness of 
our integration method regarding several aspects, at least 
for physical conditions similar to the investigated case of a 
one-dimensional CCSN: Besides convergence regarding the 
spatial resolution the test evinces that the mixed explicit- 
implicit time-integration scheme does not produce time-step 


dependent numerical artefacts. Furthermore, the test justi¬ 
fies the time-explicit treatment of the hydrodynamic quan¬ 
tities in the source-term integration. Finally, the test verifies 
that the velocity-dependent terms dropped in Eqs. |t| are 
truly insignificant. 

Two additional calculations with different closure pre¬ 
scriptions are conducted: In one run, “Mi closure”, the clo¬ 
sure in Eq. (28b I and in the other run, “maximum entropy 
closure”, the closure in Eq. (28d I is used to express the Ed¬ 
dington factor x■ For simplicity, we keep using the Minerbo 
prescription for the 3rd-moment factors q, cf. Eq. (331. Since 
the 3rd-moment terms are of minor relevance compared to 
the terms including the 2nd moments (recall that the 3rd- 
moment terms vanish in the energy-integrated lst-moment 
equation) this is a justified approximation to test the dom¬ 
inant impact of using a different closure prescription. All 
three simulations that make use of different closures give 
almost identical results, as is shown in Fig. [Tl] for selected 
quantities. This important result is consistent with the find¬ 
ings in Sec. |4.3.1| and it suggests that the AEF algorithm 
is sufficiently self-consistent when applied to CCSNe, in the 
sense that the outcome of the calculation is rather insensi¬ 
tive to the precise shape of the closure. 


5 SUMMARY 

We presented the neutrino-transport code ALCAR that 
was developed to perform multidimensional simulations of 
CCSNe and (different stages of) NS-mergers. The energy- 
dependent neutrino-transport scheme is based on the multi¬ 
group evolution (with full energy-bin coupling) of the first 
two angular moments of the specific intensity defined in the 
frame comoving with the fluid, and it takes into account the 
dominant 0(v/c) terms describing fluid-advection, aberra¬ 
tion and Doppler shift. The resulting system of equations 
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for the neutrino energy density and the three components 
of the flux density (cf. Eqs. 0 ) is closed by an approxi¬ 
mate prescription for the Eddington tensor, which assumes 
that the specific intensity is axisymmetric around the di¬ 
rection of the flux-density vector and that the remaining 
single parameter \ is given by an algebraic function of the 


evolved radiation moments (cf. Eqs. (281). Thereby, the re¬ 
sulting AEF method circumvents the computationally de¬ 
manding task to solve for the detailed angular dependence 
of the specific intensity as it is done in Boltzmann-solvers. 
In contrast to the standard FLD method, the AEF method 
consistently evolves the 1st moments and, hence, it is po¬ 
tentially more accurate (see, e.g., the test in Sec. 4.3.1|), i t 
allows radiation anisotropies to be described (cf. Sec. 4.2.11, 
and it ensures the conservation of the total (radiation plus 


fluid) momentum and energy up to 0{v/c) (Baron et al. 
1989| Cernohorsky & van den Horn 19901. Finally, a com¬ 
putationally relevant difference between the AEF and FLD 
methods is that the time step in the case of a time-explicit 
advection scheme is required to be considerably smaller in 
the FLD than in the AEF formulation, in practice forcing 
realizations of the former method to employ fully implicit 


time-integration algorithms (cf. Sec. 2.41. 


Our numerical scheme essentially follows the ideas pre¬ 
sented in Pons et al. (20001 and Audit et al. (20021, which 


have recently been implemented also in a number of photon- 
transport (e.g. HayesjkNormau] 2003jJAubert&Teyssier 


2008||Sijdowski et al.|2013||Skinner fc Qstriker|2013 McKin 


ney et al. 2014I and neutrino-transport (Shibata & Sekiguchi 


2012 Takahashi et al.|2013 Q’Connor|2014 1 codes. The ba¬ 

sic strategy is to utilize the hyperbolic nature of the moment 
equations to employ a Godunov-type scheme, in which the 
advection fluxes between grid cells are given as solutions of 
Riemann problems. Thanks to a neutrino-number conserva¬ 


tive scheme developed by Muller et al. (2010) to handle the 


Doppler shift terms we avoid the simultaneous evolution of 
the number-related moments together with the energy-based 
moments. A distinctive feature of the presented scheme com¬ 
pared to existing FLD and most Boltzmann-type neutrino- 
transport solvers is that all except the source terms are inte¬ 
grated explicitly in time. Although in this case the time step 
resulting from the Courant condition is comparable to the 
light-crossing times of single grid cells, it will usually only 
be marginally (or at least not several orders of magnitude) 
lower than the already small fluid time steps in the hot and 
dense physical systems this code is supposed to be applied 
to. To capture the transition to the stiff parabolic limit of 
the moment equations in a numerically consistent and sta¬ 
ble fashion, the source terms are handled time-implicitly and 
the upwind-type Riemann solver switches to a central-type 
solver in optically thick regions. Since the implicitly han¬ 
dled source terms are only functions of the local neutrino- 
gas properties, a computationally convenient feature of this 
explicit-implicit integration scheme is that the method can 
be parallelized with high efficiency. 

We conducted a series of tests to assess the quality of the 
AEF method and to check the correct implementation of the 
velocity-dependent terms, the source terms and the coupling 
to the hydrodynamics part of the code. By means of one- and 
two-dimensional test problems it was shown that the AEF 
method allows for a stable and self-consistent evolution of 
the radiation field in the full range between isotropic diffu¬ 


sion and free-streaming, including the accurate description 
of frame-dependent effects such as Doppler shift and diffu¬ 
sion in static and moving media. Although this was done 
here in two dimensions, the code readily generalizes to three 
dimensions. 


Two additional tests specifically focused on (one¬ 
dimensional) neutrino transport in CCSNe. In the first test 
the hydrodynamic background, consisting of a proto-NS con¬ 
figuration, was held fixed to compare the neutrino fields re¬ 
sulting from an AEF scheme with different Eddington fac¬ 
tors with the outcomes of an FLD scheme with different 
flux-limiters and of a more accurate Boltzmann scheme. 
The essential findings were that the AEF solvers repro¬ 
duced the results of the Boltzmann solver slightly more 
accurately than the FLD scheme and that using different 
closure prescriptions has less impact on the solution in an 
AEF scheme than in an FLD scheme. In the last test we 
performed a fully dynamic core-collapse simulation up to 
300 ms post bounce and we found very good agreement with 
the corresponding results obtained with the Boltzmann- 
type RHD codes VERTEX-PROMETHEUS and AGILE- 
BO LTZ TRAN. For this scenario we conducted additional 
test runs which checked the robustness with respect to our 
integration algorithm and that revealed the convenient out¬ 
come that the physical results only marginally depend on 
the actual choice of the closure relation. 


Although in this paper we investigated many cases in 
which the AEF method yields results comparable to the 
Boltzmann equation, one should keep in mind that the com¬ 
putational advantages of the AEF method compared to a 
Boltzmann solver do not come for free. That is, the closure 
relation between angular moments cannot be fulfilled to the 
same degree for arbitrarily complex radiation fields. A re¬ 
lated, particular shortcoming of the AEF method is that 
even in the optically thin limit radiation fronts interfere 
with each other, which is an immediate result of the clo¬ 
sure being a generally non-linear function of the evolved 
moments (for consequences of this features in the case of 
a post-merger BH-torus system, see the appendix of [Just] 
et al.| 20151. Nevertheless, we consider these deficiencies to 


be not overly restrictive for our purposes since the present 
code is primarily designed to describe systems in which a 
single extended source dominantly determines the evolution 
of the radiation field. 


We have already started to operate the described code 
to examine the combined neutrino- and magnetic-field ef¬ 
fects in two-dimensional CCSNe ([Obergauli nger et al.|2014| ) 
and to study the impact of neutrino transport on outflows 
from post-NS merger BH-accretion tori (Just et al. 20151. 
Since here we only discussed test problems which are sim¬ 
plified in one way or another we refer the reader to these 
mentioned papers for more specific discussions of results ob¬ 
tained with the AEF scheme in multidimensional applica¬ 
tions. In the future we plan to improve the presented code 
by supplementing the coevolution of /r- and r-neutrinos, re¬ 
fining the set of neutrino-interaction channels, and adding 
relativistic corrections. 
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APPENDIX A: DERIVATION OF THE 
3RD-MOMENT FACTOR FOR THE MINERBO 
CLOSURE 


Here we outline the derivation of Eq. <[33j), which expresses 
the 3rd-moment factor q as a function of the flux factor / 
assuming that the radiation field obeys the Minerbo closure. 
The Minerbo closure can be found by maximizing the en¬ 
tropy functional for a particle configuration with a low num¬ 


ber density (Minerbo 1978 see also Cernohorsky & Bludman 
|1994[). This leads to the axisymmetric distribution function 




(Al) 


in which fj, is the cosine of the angle between the direction of 
particle momentum and the symmetry axis, and the param¬ 
eters a, rj are two Lagrange multipliers. The latter can be 


eliminated using the definition of the 0th and 1st moments 
(cf. Eq. {5}), 


Ms)7 daF= (/rc) 3 ?”’ h “' 


(A2a) 
sinh a \ 

a J ’ 

(A2b) 


which results in 


F 1 

f = — = coth a - = L(a) 

cE a 


(A3) 


and, hence, in a = L~ 1 (f), where L and IT 1 are the 
Langevin function and its inverse, respectively. Since L can¬ 
not be inverted analytically, we employ for all practical pur- 
poses the following approximation of IT 1 Cernohorsky & 


Bludman (1994): 


a = L~ 1 (f): 


15/ 


5-3 p + p- 3/ 4 


(A4) 


of which the error is at most a few per cent. Now, the 
2nd-moment factor Xi Eq. (28aI, and 3rd-moment factor 
g, Eq. (331, are obtained from their definition in Eqs. (291 
and (311, respectively, using Eqs. (Al 1 (A41 and n- tif = /x. 


APPENDIX B: MULTIDIMENSIONAL 
CHARACTERISTIC WAVE SPEEDS 

In this section we provide supplementary information about 


the speeds A±,o (cf. Sec. 3.4 1 of the characteristic waves as¬ 


sociated with the hyperbolic part of the multidimensional 
two-moment system closed by an algebraic Eddington ten¬ 
sor. Without loss of generality, we consider characteristic 
waves propagating along the x-direction and assume that 
the lst-moment vector F= (F x ,F y ) = (y, \J\ — U 2 )\F\ lies 
in the x — y plane and forms an angle qf, where /r = cos «f, 
with the x-axis. The wave speeds are then obtained as eigen¬ 
values of the Jacobi matrix whose components are given in 
terms of a general closure x(e, f ) as follows: 


dF x dF x dF x 
dE ’ dF x ’ dF v 


= (0,1,0) 


(Bla) 


2 dP xx _ c 2 
C dE ~ 2 

2 t)P xx _ c/r 

c ~dFj - 27 

2 t)P xx _ c^l U 


l- M ‘ + (3p - 1 ) \ e fo~ f df +x 


dx ,-°X 


/fy ( 3 m 2 -1)-2(/x 2 -1)(3 X -1) 


dF v 


2 / 


/|y(3^-l)-2/r 2 (3.X-l) 


2 dP X y _ (?U\J 1 - M 2 
c dE 2 

2 t)P X y _ C^l [P 




dF x 


2 / 


3 X -1 + M 2 (2 + 3/0-6 X 


(Bib) 

1 

(Blc) 

(Bid) 

(Ble) 

)]• 

(Blf) 
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Figure Bl. Comparison of characteristic speeds A± s o/) (with A— < Ao < A+) as functions of the direction cosine /x (left) and flux 
factor / (right) for the Minerbo, M\ and Monte-Carlo closure (cf. XMinerbcnXMi and XJanka in Eqs. ( |28| >, respectively). The thick lines 
denote the exact wave speeds while the thin black lines depict for the Minerbo closure the linear expansion in /x, cf. Eq. |4b|. 


lOPxy 

dF v 


c/i 

V 


3/^(1 -/x 2 ) + (2/x 2 -1)(3 X -1) 


(Big) 


We plot in Fig. |Bl| the resulting wave speeds as functions of 
/x and / for three different closures together with the linear 


expansion in /x, Eq. (461, for the Minerbo closure. It can be 


seen that all considered closures lead to rather similar wave 
speeds. Moreover, the /x-dependence of the wave speeds is 
very close to linear, which suggests that using the linearized 
versions instead of the exact wave speeds as signal speeds for 


the Riemann solver, cf. Eq. (481, is a justified approximation. 


APPENDIX C: NUMERICAL METHODS USED 
FOR THE COMPARISON IN SEC. 4.3.1 

Cl Flux-limited diffusion solver 

To construct the spherically symmetric FLD solver used to 
find the corresponding solutions in Sec. |4.3.1| we start from 


the presented AEF scheme and essentially drop the evolu¬ 
tion equation for the radial flux density F r = F and all 
velocity-dependent terms. What remains to be determined 
is a suitable numerical representation of the flux Ffld (cf. 
Eq. ( |22| )) at each cell interface. To this end, we first com¬ 
pute the cell-centered version of the flux, for which we need a 
cell-centered representation of d r E. The latter is discretized 
as 

d r E —> (E\+i - E i _ 1 )/(2Ar i ), (Cl) 

where i denotes the radial grid index and An = r -, i — r- i. 

'22 

Out of the resulting cell-centered values Ffld , i of the flux 
density, we compute the cell-interface values as 

^FLD.i+i = / ^FLD,i+i^' FLD . i + (l~^FLD,i+i)^I i 'LD,i+l , (C2) 

where the interpolation parameter Afld £ [0,1] is intro¬ 
duced to ensure that the numerical method is based on 
central-type fluxes in the parabolic regime, where Afld —> 
1/2 should hold, and on upwind-type fluxes in the hyper¬ 
bolic regime, where Afld —► 1 for / = Ffld/(c.E) —» 1 
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and Afld —> 0 for / —> — 1 should hold. We compute 
Afld as a function of the signed, average flux-factor f i+ i = 

{Fflv;JE\ + Fffy>,\+i / E\ +1 ) / {2c) as 


A 


FLD,i+ 


i = max{ min{ / i+ i + 1/2,1}, 0} . 


(C3) 


The time step At for the explicit integration is chosen to 
fulfill At < minj{rFLD.i}, cf. Eq. (26b I with a ss 0.5. 


C2 Tangent-ray Boltzmann-solver 


For calculating the reference solution in Sec. |4.3.l1 we employ 
a time-independent tangent-ray variable-Eddington-factor 
(TR-VEF) scheme closely analog to what is described in 
Chap. 83 of |Mihalas &; Mihalas (1984 1 and in Rampp & 
Janka ( 2002| ). For the discretization of the two-moment sys¬ 
tem (cf. Eqs. 0 with dt = 0, v = 0 and the source terms 
expressed as in Eq. (141) we use a finite-difference scheme in 
which we interprete the energy densities E, to be located at 

cell centers and the flux densities F- , i to be located at cell 

' 2 

interfaces. Using first-order differences for the radial deriva¬ 
tives, this leads to the following linear system of equations: 


Pi = n (i = l,...,JVr). We employ a finite-difference pre¬ 
scription in which we interprete j\ to be located at cell cen¬ 
ters and h, , i to be located at cell interfaces of the corre- 
2 

sponding tangent-ray grid. The linear system of difference 
equations to be solved then follows in an analog fashion as 


for the two-moment system, cf. Eq. (C4 1 . For the numeri¬ 


cal angular integration of the specific intensity we use the 


quadrature weights of Yorke (1980). 


In practice, we obtain the initial values of x by solving 
the Boltzmann equation with E = F = 0. Subsequently, 
the two-moment system and the Boltzmann equation are 
iteratively solved until convergence is achieved. 


F-. i — F. 


2 F, 


Xi+iEj+i — Xi-Ej 

Ar,, i 


+ CK a ,i-Ei = C/t a ,iE i eq , 

E a 


(C4a) 


-C'i+l / \ 

4-— (3x i+ i — l) 


+ 4 “ 1 ) +«tra,i+I^+i =0 


(C4b) 


for i = 1,..., N r , where Ar., i = n+i — n and E -, i, F\ are 

' 2 ' 2 

linear interpolations of the nearest neighbors on the corre¬ 
sponding staggered grid. The Eddington factors x required 
to solve Eqs. ( |C4[ ) have to be obtained from the Boltzmann 
equation. For the considered type of interactions, the latter 
reads (e.g. Cernohorsky et ah|l989 1 


lid r T{r, n) + 


1-/U 


d^X{r, p) = 


KaT eq + — (kS c E + kI m F) - (Ka + «5)X(r, fi) , (C5) 

where p is the cosine of the angle between the radial di¬ 
rection and the direction of neutrinos in momentum space, 
X eq is the specific intensity corresponding to neutrinos 
being in thermal equilibrium, and the scattering opacity 
k s = k° — K q /3 is decomposed into an isotropic (super¬ 
script “0”) and an anisotropic (superscript “1”) contribu¬ 


tion. For solving Eq. (C51 we first make a change of variables 
(r, p) —> (s = fi r, p = x/1 — p 2 r) and we use 

/( S ,p) = (I(p)+T(-p))/2, (C6a) 

h{s,p) = {X{^-X {-^))/2 (C6b) 


to rewrite Eq. ( |C5[ ) into the following system of equations 

d s h(s,p) + (Ka + K°)j(s,p) = K a I eq + K°J- E , (C7a) 

47T 


dsj(s,p) + (K a + K°)h{s,p) = K S J— pF. 

47T 


(C7b) 


Equations ( |C7[ ) are independently solved along each 
tangent-ray characterized by a constant impact parameter 
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